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Abstract. We study the effective mass of the bipolarons and essentially the possibility to get both light 
and strongly bound bipolarons in the Holstein-Hubbard model and some variations in the vicinity of the 
adiabatic limit. Several approaches to investigate the quantum mobility of polarons and bipolarons are 
proposed for this model. First, the quantum fluctuations are treated as perturbations of the mean- field 
(or adiabatic) approximation of the electron-phonon coupling in order to calculate the bipolaron bands, ft 
is found that the bipolaron mass generally remains very large except in the vicinity of the triple point of 
the phase diagram (see |2^), where the bipolarons have several degenerate configurations at the adiabatic 
limit (single site (SO), two sites (SI) and quadrisinglet (QS)), while the polarons are much lighter. This 
degeneracy reduces the bipolaron mass significantly. Next we improve this result by variational methods 
(modified Toyozawa Exponential Ansatz or TEA) valid for larger quantum perturbations away from the 
adiabatic limit. We first test this new method for the single polaron. We find that the triple point of the 
phase diagram is washed out by the lattice quantum fluctuations which thus suppress the light bipolarons. 
Further improvements of the method by hybridization of several TEA states do not change this conclusion. 
Next we show that some model variations , for example a phonon dispersion may increase the stability of 
the (QS) bipolaron against the quantum lattice fluctuations. We show that the triple point of the phase 
diagram may be stable to quantum lattice fluctuations and a very sharp mass reduction may occur, leading 
to bipolaron masses of the order of 100 bare electronic mass for realistic parameters. Thus we argue that 
such very light bipolarons could condense as a superconducting state at relatively high temperature when 
their interactions are not too large, that is, their density is small enough. This effect might be relevant for 
understanding the origin of the high Tc superconductivity of doped cuprates far enough from half filling. 

PACS. 71.10.Fd Lattice fermion models (Hubbard model, etc.) - 71.38.-|-i Polarons and electron phonon 
interactions - 74.20.Mn Nonconventional mechanisms (spin fluctuations, polarons and bipolarons, res- 
onating valence bond model, anyon mechanism, marginal Fermi liquid, Luttinger liquid, etc.) - 74.25.Jb 
Electronic structure 

1 Introduction are described by a hard core boson model which could 

become superfluid. Unfortunately, their calculations also 
1.1 Specific Problem for high TC Superconductivity show that when the electron phonon couphng increases 
$H beyond the Migdal instability, the effective mass of these 

- ■ Superconducting materials at temperatures significantly bosons grows exponentially fast and becomes so huge that 
higher than the maximum Tc predicted by MacMiUan f hopeless to get superconductivity m this model, 

r ; J J T->/-ic J i.- -J. 4-- 1 at least at non-negligible temperatures, 

for the standard BLb superconductivity are exceptional. ° ° ^ 

Up to now, there is a wide variety of such materials which One of us(SA) conjectured in that the interplay 
are all cuprates built with Cu02 planes and with many of the electron-phonon coupling witn a direct electron- 
kinds of interlayer dopants ||] . electron repulsion could reduce significantly the bipolaron 
When the electron-phonon coupling increases too much, effective mass and thus favor high Tc superconductivity, 
it is known for several decades (Migdal 0) that the BCS In the present paper we support this conjecture by cal- 
theory should break down because of lattice instabilities, culating the effective mass of a single bipolaron in a 2d 
These instabilities are associated with the formation of model which involves both an electron-phonon coupling 
polarons and bipolarons. Alexandrov et al. developed and a direct electron-electron repulsion. For this purpose 
later a theory of bipolaronic superconductivity where the we choose to investigate first this effect in the Holstein- 
electrons form strongly bound on-site bipolarons. They Hubbard model because of its simplicity. 
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In the absence of Hubbard repulsion, we confirm that 
the effective mass of the bipolaron is indeed very large 
Ip^ , which is incompatible with a high Tc superconduct- 
ing phase. When the Hubbard term is increased new bipo- 
laronic states become stable. They are 2-site bipolarons 
which consist of two neighboring polarons bound by their 
magnetic interaction in a singlet state and also a bipo- 
laron called a "quadrisinglet" (QS) which consists of the 
combination of four singlets sharing one common site. In 
the parameter region where these bipolaronic states have 
nearly degenerate energies, the effective bipolaron mass 
is sharply reduced. The drastic mass reduction is due to 
resonance between the different bipolarons. Certain vari- 
ations of the model, such as a phonon dispersion, might 
even increase the binding energy of the bipolaron while 
allowing it to keep a very light effective mass for realistic 
parameters. 

We first discuss some early known results about the 
effective masses of polarons and bipolarons. Adiabatic re- 
sults (section 1^) described in ||2^ are briefly recalled. Next 
we treat the quantum lattice fluctuations as a perturba- 
tion of the adiabatic limit. The main effect is to lift the 
bipolaron degeneracy both due to translational invariance 
and to the possible existence of several kinds of adiabatic 
bipolarons with almost the same energy (section ^ . This 
correction is only valid for a very small quantum lattice 
parameter. 

To extend the field of validity of our calculations we 
next propose to use the Toyozawa variational form. Quan- 
tum polarons and bipolarons are approximated by a self- 
consistent Bloch wave that is exact at the adiabatic limit. 
At this limit it is demonstrated in |Q that there is not a 
great loss of accuracy if the shapes of polaron and bipo- 
laron are exponentials. Thanks to this approximation we 
gain much simplicity for the variational form. We first ap- 
ply this method to polarons (section U) and next to bipo- 
larons (section ^ . For very small quantum lattice fluctua- 
tions the results obtained by perturbation of the adiabatic 
limit are practically recorvered but there are significant 
deviations when the fluctuations become larger. First or- 
der transitions that cannot exist physically are washed out 
by hybridizing several Toyozawa variational forms. Actu- 
ally, the regime where the bipolaron mass is sharply re- 
duced is swallowed up in the domain where the ground 
state is unbound polarons. We demonstrate in section (||) 
that this undesirable phenomena can be avoided by varia- 
tions ofthe model which increase the stability of the (QS) 
bipolaron (e.g. a phonon dispersion with the appropriate 
sign). 



where j and k represent lattice sites, T is the trans- 
fer integral between nearest neighbor sites < j,k >. The 
electrons are Fermions represented by the standard anti- 
commuting operators Cj'^ and Cj,a at site j with spin 

(T =t or J,. Hi^a- = C^^^Cj^cr and = -t- n^j^. a'^^ and aj 
are standard creation and annihilation boson operators of 
phonons and hujo is the phonon energy of a dispersion-less 
optical phonon branch, g is the on-site electron phonon 
coupling and v the on-site electron-electron repulsion (Hub- 
bard interaction). 

We choose Eq — Sg^/Ticjo as energy unit as in [^ . 
Defining the position and momentum operators as 



hujQ 



hujo ^ 



with the commutation relation [uj^pj] 
we obtain the dimensionless Hamiltonian 



(2) 

(3) 

(4) 



^^2 



(5) 



<j,k>,cr 



Our reduced dimensionless parameters are 
T „ 1 

t ■ 



Eo = Sg'/hLOo 

T 
E^ 
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(6) 



Despite the primitive nature of our model, it may catch 
important aspects of reality. However, we shall also demon- 
strate at the end of this paper that certain model varia- 
tions could be favorable for bipolaron mass reduction. For 
example, we may introduce a coupling between nearest 
neighbor atoms so that the phonon branch is no longer 
dispersionless. Then the new Hamiltonian is the sum of (|g) 
and the extra energy term —C J2<i j> UiUj where < i,j > 
represents all the pairs of nearest neighbor sites i and 
J. When C > 0, the bipolaron mass reduction is en- 
hanced while it remains strongly bound. Thus we demon- 
strate that relatively minor changes in the model may fa- 
vor (or disfavor) superconductivity at relatively high tem- 
perature. 



1.3 Polaron and Bipolaron Effective Mass 



1.2 The Holstein-Hubbard Model 

Let us first recall our notations for the model we study 
here. Its Hamiltonian is 

+gnj (a+ + a^) + vrij^^Uj^i (1) 



Let us first briefly recall some standard results about the 
effective masses of polarons and bipolarons. The Lang- 
Firsov unitary transformation yields a new Hamilto- 
nian 
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where Hlf = e-"-^^" He''^'^'' with Slf = \ Y.jPj'^j 

After this transformation, the creation operator Cj*"^ 
at site j acts on the vacuum by creating both an electron 
and a lattice distortion 



+ p-»Pj/2| 



> 



(8) 



that is, a polaron. A standard but rough mean-field ap- 
proximation consists in taking the average of the transfer 
integral for unperturbed phonons. We obtain an approxi- 
mate formula for the transfer integral of this polaron 



masses are so huge that it is unreasonable to consider that 
they could exhibit a Bose condensation |p^ . 

We perform here a numerical calculation of the effec- 
tive mass of the bipolaron (and also the polaron) in order 
to show that in some specific regions of the parameter 
space, when the Hubbard term becomes comparable with 
the electron-phonon binding energy, these effective masses 
can be drastically reduced so that Bose condensations of 
bipolarons become plausible. 



(9) 2 The Mean-field Holstein-Hubbard Model 



When there is a single electron in the system the electron- 
electron interaction does not play any role. Then the ef- 
fective mass of a single polaron is defined as the inverse 
of the second order derivative versus the wave-vector q 
of the polaronic energy E{q), that is, Tlf- The effective 
mass of the polaron is that of the bare electron multiplied 
by exp [l/(8a)]. It becomes exponentially large when a is 
small. 

This approximation tends to become right only when 
the operator corresponding to the transfer integral in eq.^ 
has small fluctuations. This condition is fulfilled when the 
pre-factor t is small. It is also fulfilled when a is large, 
that is, for a weak electron-phonon coupling g compared 
to the phonon energy hiVQ. In the antiadiabatic limit {a 
large), the model becomes a Hubbard model with an on- 
site effective electron-electron interaction U = (C/ — 1/4) 
which is attractive for t7 < 1/4 and repulsive for [/ > 1/4 
and where the transfer integral has been renormalized. 
The negative U model is expected to have superconduct- 
ing phases [|l4[ for non- vanishing band filling. However, we 
treat here the opposite case a small which is close to the 
adiabatic limit. 

Actually when t is small, our numerical results agree 
with formula (^V For larger t, the effective mass of the po- 
laron given by (|9|) becomes larger than the mass we com- 
pute. Note that our result should be more reliable because 
it yields a lower variational energy for a single polaron. 

The effective mass of the bipolaron has been calcu- 
lated by Alexandrov et al in the same limit {t small) [0 
for strongly bound bipolarons (that is, for U small) and 
considering the kinetic energy term in eq.0 as a pertur- 
bation. In our dimensionless units they found the transfer 
integral th for a bipolaron 



tb 



At' 



.-l/(2a) 



1-4C/ 



(10) 



If one extrapolates naively this formula for larger U , 
one would find that % becomes infinite. Of course, the 
associated effective mass of the bipolaron cannot vanish, 
but our results nevertheless demonstrate that it is sharply 
depressed not far from the region U « 0.25. Comparison 
of formula (^) and ^ shows that in the region where both 
a and t are small the effective mass, of the bipolaron is 
much larger by many order of magnitude than the polaron 
mass which is itself much larger than that of the bare elec- 
tron. In most physically realistic situations, the bipolaron 



We calculate first the adiabatic bipolarons which are ground- 
state of a mean-field Hamiltonian. They are the exact so- 
lutions in the adiabatic limit when 7 = is zero (that 
is, when the atomic kinetic energy is negligible). These 
spatially localized solutions are degenerate under lattice 
translation. This degeneracy is lifted when the atomic ki- 
netic energy is taken into account. Within a perturbative 
treatment, this explicitly gives bands of extended bipo- 
larons characterized by a wave-vector. The inverse cur- 
vature of the lowest band at zero wave-vector yields the 
bipolaron effective mass. This calculation become exact in 
principle in the limit of small 7. Note that similar meth- 
ods were already developed in |22] to calculate the effective 
masses of discommensurations in Charge Density Waves. 



2.1 The adiabatic regime 

The Hamiltonian eq.^ can be written as the sum of three 
terms 

H = H,i+Hph + Hf (11) 

where Hei and Hph are decoupled electron and phonon 
Hamiltonians respectively and Hf is a fiuctuation term. 



<i,j>,a- 



i 



(12) 
(13) 

(14) 



fii and Ui are variational parameters which are deter- 
mined by minimizing the ground-state energy of the ef- 
fective Hamiltonian Had — H^i + Hph- It comes out that 
Ui — (ui) and fii — (rii) are the average of the correspond- 
ing operators. The standard mean-field approximation for 
the electron phonon coupling consists in neglecting the 
fiuctuation energy Hf. 

Minimizing the ground-state energy of Hamiltonian 
eqjl^ also yields 

(u,) - -(n,)/2 (15) 
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Then, the ground-state of the mean field Hamiltonian 
Had has the form 



> (16) 



A pair of electrons with the electronic wave function 
ipij is created as well as the corresponding lattice distor- 
tion Ui. The electronic wave function is a singlet state, 
that is, a symmetric function of (i, j): = ipj.i ■ K ful- 
fills an extended nonlinear Schroedinger equation which is 
exactly the same as in the adiabatic case at a = p5[. 



t 



+ UStj Tp.,J = Eel^puj (17) 



A is the four-dimensional discrete Laplacian and Ui — 
""■1,1 + "-ij with 

n^A^y,H^.J\' (18) 



j 

The square root of the mean square lattice fluctuation 
{{ui — Ui)'^)^^'^ = a = 7^/^ is small of order a. Thus, the 
mean-field approximation obviously becomes exact in the 
adiabatic limit a — > when there are no lattice fluctua- 
tions. 



2.2 Bipolarons from Anti-integrable limit and 
variational approximations 

For an easy understanding the reader should refer to our 
early paper were the adiabatic (or mean field) bipo- 
larons were investigated in detail in the two-dimensional 
model by continuation from the anti-integrable limit {t = 

0)- . _ . 

The main result of pq | is that we found a quite rich 
phase diagram with first order transition lines in the pa- 
rameter space (U,t). For large t the electrons remain ex- 
tended and do not self localize as bipolarons. For small 
t there are several kinds of structures that compete to 
be the bipolaron ground-state. These bipolarons were de- 
noted (SO), (SI) and (QS). Bipolaron (SO) is mostly local- 
ized at a single site and has square symmetry. Bipolaron 
(SI) consists into a bound pair of polarons in a magnetic 
singlet state localized on two neighboring sites. It breaks 
the square symmetry and is oriented either in the x di- 
rection (Sl)x or the y direction (Sl)j,. The quadrisinglet 
bipolaron (QS) is a combination of four singlet states with 
a common central site and has square symmetry. 

Interesting properties are obtained at a triple point 
corresponding to the intersection of three first-order tran- 
sition lines. At that point, and apart from the transla- 
tional degeneracy, there are four different degenerate bipo- 
larons (SO), (QS), (Sl)^ and (Sl)^. We shall see that the 
quantum lattice perturbations hybridize these degenerate 
states and hence drastically enhance the bandwidth of 
the bipolaron or, equivalently, reduce its effective mass. 
Within a classical picture we already noticed that the en- 
ergy barrier ( Peierls-Nabarro barrier) which has to be 



overcome to move the bipolaron through the lattice was 
drastically reduced. 

We also investigated in p5| some approximations with 
exponential variational forms for the bipolarons that al- 
low analytical calculations. The exact phase diagram cal- 
culated numerically was reproduced with a quite good ac- 
curacy with the following forms 



(19) 



ix-l| + |ial + bx| + |jj,|) 



(20) 



c 



(bx| + |i„|)/^(|^.±l| + |i„|) 



|ix| + |iy±l|: 



'(A- 



^;^(b.l + b„±l|)) (21) 



for bipolarons (SO), (SI) and (QS) respectively. A, B, and 
C are normalization factors and the parameters A, Ai and 
A2 are optimized for energy minimization. We shall de- 
velop here a quantum analogous version of these approxi- 
mations to improve our methods. 



3 Quantum Lattice Corrections 

We now treat the mean- field fluctuation Hf = 1/2 X]i("i^ 
Ui){ni — fii) as a perturbation that lifts the translational 
degeneracy of the mean field bipolarons ( [l^ ) , whose wave 
functions are denoted \Q^{j) > where S represents bipo- 
larons (SO), (Sl)a;, (Sl)j^ or (QS). The index j is the site 
where the bipolaron (S) is located. For bipolarons (Sl)a; 
and {Sl)y which occupy two adjacent sites {jx,jy) and 
Ux + l,jy) or {jx,jy) and {jx,jy + 1), respectively, we 
choose by convention j = {jx,jy)- To treat the mean- field 
fluctuation in lowest order, the initial Hamiltonian ( |Tl| ) 
should be projected and diagonalized in the subspace gen- 
erated by all these translated wave functions. We already 
noticed that the bipolaron energies might be degenerate 
or almost degenerate so that we should take into account 
their possible hybridization. The eigenstates should have 
the general form 

l^>=E«^^^l^''0')> (22) 

where as,j are coefficients to be determined by extremal- 
ization of < f2\H\f2 > with the normalization constraint 
< n\n >= 1. Both < n\{2 > and < i7|H|i7 > are 
quadratic functions of asj 

< [2\H\[2 > = E a*s,^M(^s\^),iS',J)as'.J (23) 
<n\Q>= J2 a*s^,P(^sMS',j)as',j (24) 

where the coefficients of matrices P and M are defined as 

P{s,^).iS'.J) = < f2s{^)\ns'U) > (25) 
M(s,,),(S'.j) = < ns{t)\H\ns'{j) > (26) 
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It is important to take into account the fact that the eigen- 
states of the self- consistent mean-field Hamiltonian Had 
are not orthogonal one with each other, since the matrix of 
scalar products is not diagonal. For two normalized wave- 
functions |J7 > and |J7' > with the form ( p^ and with 
electronic wave functions {ipij} and j}, electronic den- 
sities fii = ~2wi and h'^ — —2u[ respectively the scalar 
products eqs.(^6|),(^) can be calculated explicitly for the 
Hamiltonian IM). 



< n\n' > 




N- < > 



(27) 



exp - 



- fe 



n i,j 



4^ 



+U J2 V^mV^M - 2 E i^lj^^j)] (28) 

i i,3 

where ^V'ij' = V'fe,; is the discrete Laplacian on a 4d 
lattice. 

The extremalization equation of < [2\H\fi > with re- 
spect to A = {as,j} with the normalization condition 
A*. P. A = 1, is M.A = E.P.A where E is the Lagrange 
parameter, which is also the eigenenergy. It can be writ- 
ten as an eigenvalue problem for the normalized vector 
B = P^/^A 



P^/^.M.P^/^B = E.B 



(29) 



Note that the extensive term < f2\f2' > in the 
second term of eq. (p8|) does not disturb the calculations. 
It yields a constant term in the effective Hamilto- 
nian P~^^^.M.P~^^^ which is nothing but the zero point 
phonon energy of the system with size N (without elec- 
trons). 

Because of the translation invariance of the model. 



M, 



) and P(s.i).(fi'.-i) only depends on j 



iS.i),(S',j) 

As a result, eq.^ can be diagonahzed as combinations 
of plane waves with the form Bsj{K) — J2s ^siK) 
with wave vector K which fulfills the eigenequation 

P"^/2(if).M(if).P"^/2(if).B = E^iK).BiK) (30) 
with the Fourier coefficients 



(Sj%(S'j+«)e 



iKr, 



Kr. 



(31) 
(32) 



Then the diagonalization of the 4x4 matrix 
'P~'^/'^{K)M.{K).'P~'^/'^{K) yields the eigenenergies 
Ei,{K). Figure (|l|) shows an example of calculation of 
these four bands in the vicinity of the triple point. Thus 
when there are four bipolarons that are metastable (e.g. 
in the vicinity of the triple point of the phase diagram) 
one obtains four bipolaron bands E^{K). Within our ap- 
proach the number of bipolaron bands is equal to the 
number of metastable states for the adiabatic bipolaron 
which provides the base about which we expand the eigen- 
states. In other regions of the phase diagram the number 
of metastable bipolarons changes, which induces (unrealis- 
tic) discontinuities for the number of bands. For example, 
when U — only the bipolaron (SO) is metastable, and 
there is only one bipolaron band. 

However the lowest bipolaron band does not exhibit 
very sharp changes despite a small discontinuous varia- 
tion. The reason that the upper bands are not reliable is 
that the energies of these states might be also degener- 
ate with phonon excitations of the bipolaronic states of 
the lower band. The real excited states involve complex 
hybridization between these states. 

Conversely, the bipolaronic states with the lowest en- 
ergies should not hybridize significantly with the higher 
energy states involving phonon excitations. Thus we con- 
sider that the lowest-energy bipolaron band provides a 
reliable description of the bipolaron excitations close to 
its ground state. We use it to measure the bipolaron ef- 
fective mass, that is, the inverse of the curvature Tf, at 
zero wave vector = 0. Tf, is constant in all directions 
because of the square symmetry (figj^). It can be viewed 
as the effective hopping coefficient for the bipolaron tun- 
nelling through the lattice and can be compared with the 
prediction of jl^ given by formula (^0|), which is valid at 
both U and t small. 

Fig.H exhibits the ratio Tf, /if, as a function of the ef- 
fective transfer integral t for U = for several values of 
the quantum parameter a. For U small this ratio goes to 
1 when t goes to 0, which confirms the validity of formula 
( |l0|) predicted by Alexandrov et al. |17 in that regime. 
We also note that beyond this regime when the parame- 
ters {U, t) are larger than 0, Ti, becomes significantly larger 
than ib, or equivalently the bipolaron effective mass cal- 
culated numerically drops faster than predicted by (^^. 

The insert of fig.^ shows the bipolaron energy gain 
compared to a pair of free electrons |^ at L/ = and the 



fig.y shows the effective transfer integral compared to the 
bare electronic transfer integral, which is negligible at the 
scale of the electronic energy. The bipolaron effective mass 
appears to be much beyond than 10^" electronic masses 
even when its binding energy vanishes. It is clear that 
this regime C/ = is not favorable at all for the Bose 
condensation of such bipolarons that should occur below a 
critical temperature inversely proportional to the effective 
mass of the quasi-particle. 



^ The energy gain compared to a pair of free electrons is 
not an accurate binding energy for the bipolaron. The bipo- 
laron binding energy is precisely measured with respect to an 
unbound polaron state that is defined in the next section. 
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Fig. 2. Ratio Tb/tb versus t of the transfer integral numerically 
calculated and analytically predicted by formula ((7 = 

and a = 0.022 (dot-dashed line), a = 0.017 (dashed line), 
a = 0.01 (full line). 



When the Hubbard term increases for relatively small 
i, fig.§ {t = 0.04) shows that a sharp discontinuity occurs 
when the ground-state bipolaron becomes (SI) U > 0.25. 
There is a sharp increase of the tunnelling energy by five 
orders of magnitude for this bipolaron at a = 0.01. In that 
case the Peierls-Nabarro barrier calculated in the previ- 
ous paper is still very high and consequently there 
is almost no hybridization between (SO) and (SI). The 
smoothing of the discontinuity of the tunnelling energy is 
thus hardly visible. 

When t is larger, the bipolaron (QS) becomes stable for 
U w 0.23 and hybridization between (SO), (SI) and (QS) 
becomes significant. Actually the most important contri- 
bution to the tunnelling energy of the bipolaron comes 
from the hybridization between (QS) and (SI). It is re- 
sponsible for the sharp increase of the tunnelling energy 
or equivalently the sharp drop of the bipolaron effective 
mass. This quantum mobility is favored when (QS) and 
(SI) are degenerate in energy and separate by a weak 
Peierls-Nabarro barrier. Then (QS) may tunnel to one of 
the four neighboring (SI) and the latter tunnels to its 
neighboring (QS) that corresponds to the initial one trans- 
lated by one lattice spacing in the direction of (SI) and so 
on. The bipolaron tunnelling energy could reach 10~'^ the 
bare electronic energy which is not negligible anymore. 

Let us point out that such a high mobility cannot be 
obtained within the approximations used in fl^ which 
do not consider the possible degeneracies of several bipo- 
larons. The conclusions of jl^ about the physical impos- 
sibility of bipolaronic superconductivity are irrelevant for 
that situation. 



Fig. 1. Four bipolaron bands E,^{K) versus wave- vector K 
computed close to the triple point where bipolarons (SO), (SI) 
and (QS) are degenerate in energy (a = 0.017, U — 0.232, t — 
0.08). Energies increase from bottom to top. 
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Fig. 3. Ratio Tt/t versus U &t t = 0.04 and a = 0.01; insert: 
Bipolaron energy gain compared to a pair of free electrons 
versus t at U = 0. 



0.27 
U 



Fig. 4. Ratio Tt/t versus f/ at t = 0.04 (top), Tb/t versus U 
at t = 0.08 (bottom) (a = 0.01). 



4 Variational Calculation of Quantum 
Polarons 

In principle the above approach is valid only for very small 
a: that is, when the quantum lattice fluctuations are small. 
However these fluctuations may increase drastically, espe- 
cially close to the first-order transitions when there are 
several degenerate bipolarons that we are especially inter- 
ested in. Thus it is worthwhile to improve our previous 
calculations by a variational approach which should be 
equivalent to the mean field perturbation for small quan- 
tum lattice fluctuations. 



Our purpose is now to develop a quantum version of 
the variational forms used and tested in the adiabatic 
case but which could hold for larger values of a. Our ap- 
proach is a simplified version of those of Toyozawa (see 
[ p4| and [ p5) for a recent application to the polaron in 
ID). We first test this method for the single polaron and 
will extend it in the next section for the bipolarons of the 
Holstein-Hubbard model. 

Because of the invariance of the system under trans- 
lations the wave function of a quantum single polaron is 
written as a Bloch wave: 



\nP{K) >-- 



1 



|<f^(n) > 



(33) 



where yl is a normalization factor and \^^{n) > is ob- 
tained from a unique wave function |!f'^(0) > changing 
all the indices i of its electronic and atomic variables into 
i + n. This transformation is nothing but a shift of the 
wave function from site to n. 



4.1 Toyozawa approximation 

A simple variational approximation proposed by Toyozawa 
is to assume that the local wave function is similar to the 
mean-field polaron: 

E i^k-,C+).e^p (^zY,v,^,.p}j 10 > (34) 

To simplify the spin of the electron is omitted. For each 
wave vector K the variational energy 



(35) 



is a function of the scalar products which does not depend 
on j 

< f^(j)|<f^(i+p) > = exp-^^(w,+p-'y,)^ 



^ (^E^^+p^y 

i 

71 i 



(36) 



(37) 



and has to be extremalized with respect to the 2N param- 
eters corresponding to the electronic wave function {V'J'} 
and the lattice distortion {vi\. This form becomes exact 
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in the adiabatic limit and should improve the previous 
perturbation theory as it is self-consistent. 

A relation between the electronic density and the av- 
erage of the atomic displacement can be easily taken into 
account in this variational form. First, let us recall that the 
true eigenfunctions of iJ are extrema of < > in the 

full space of normalized functions ^ . For a given normal- 
ized eigenfunction f2 of H we can consider the one parame- 
ter family of normalized functions ^{S) — exp (iS.pj)\f2 > 
where the coordinate Uj of the atom j is changed into 
Uj + d. The variational energy of this wave function is 
< 1r{S)\H\1r{S) > which is equal to < 1r{0)\Hs\^{0) > 
where Hs = exp —(ipj6)H enp (ipjS) is simply obtained 
from H by changing Uj into Uj — 5. The variational energy 

< ^{5)\H\^{5) >=< [2\H\[2 > 
- i(5 (2 < n\uj\n > + < {2\nj\n >) + ^6^ (38) 

should be extremal for 6 = 0, which implies 

< n\uj\n >= < n\nj\n > . (39) 

This result is nothing but an extension to the non-adiabatic 
case of the standard relation between the average atomic 
positions and the electronic densities. 

If we now consider an extremum of < f2\H\f2 > for 
I n > normalized in the variational space defined by eq. ( |3^ ) 
and (|3^) this space is no longer invariant under the unitary 
operator e'xjp{iSpj}, but it still remains globally invariant 
under operator exp{i(5 pj } which performs a uniform 
displacement by 5 of all the atoms. We apply the same 
argument as above that is, study < n{S)\H\n{S) > where 
f2{5) = exp{i5J2jPj}^ is extremal for S = 0. This con- 
dition yields J2j < ^Wjl^^ >= -1/2 < n\nj\0 >. 
For the variational extrema with the Toyozawa form ( ^3|) 
and m, we find < '?'(0I E„ W"|>f (m) >= (E„««) < 
'P{l)\<Pjrn) >, which readily implies J^j < > = 

X^n'*^"- ^'^^ polaron, that is, for a system with only 
one electron the extremum of the Toyozawa form ( p3| ) and 
( ^ ) necessarily fulfills 

E^" = -^ (40) 

n 



4.2 Toyozawa Exponential Approximation: TEA 

Minimizing the variational form (^) for the whole set of 
2N — 1 parameters {ipi} and {vi} with condition ( pO| ) is a 
complex numerical task which moreover will become even 
more complex when extended to the bipolaron problem. 
However, we can expect that the behavior of the varia- 
tional parameters {(j^n} and will not be far from ex- 
ponential at infinity. Thus assuming simple exponentials 
for {0„} and Jw„} should not be a bad approximation as 
proposed in |2q| at the adiabatic limit. Taking into ac- 
count the normalization and condition (Eoh we postulate 



that the electronic wave function and the atomic modula- 
tion have the form: 

V^f = AAI^-I+I*"! = (1 + A2)/(1- A^) (41) 

v[ = -Bf,\^^\+^^y\ B-' = 2{^ + ^l)/{l-^,) (42) 

for each wave vector K there are only two variational 
parameters X{K) and ii{K) instead of 27V for the original 
Toyozawa ansatz which allows much simpler calculations 
although we still need a numerical minimization of (p5|). 
To that aim we use a simplex method pTf , which is the 
most efficient algorithm we tested because it avoids any 
precision problem due to the numerical computation of 
the derivatives. 

The Toyozawa Exponential Ansatz (TEA) turns out 
to be almost as good as the full ansatz when the polarons 
are small, since in that case the exponential approximates 
quite well its shape. When the size of the polaron becomes 
larger, the TEA (as well as the original Toyozawa ansatz) 
yields a first order transition. This first order transition 
is well-known to exist at the adiabatic limit a = where 
the ground state of a single electron undergoes a first order 
transition from a small polaron to a free electron [^,|6j at 
< = « 0.07486. 

We define the binding energy of the quantum polaron 
as the difference between the energy of the extended elec- 
tron at zero wave vector K = and that of the bottom of 
the polaron band. Fig|| shows the variation of the bind- 
ing energy versus t for the quantum polaron calculated in 
several different approximations including the assumption 
that: 

1. The polaron band is calculated as for the bipolaron 
bands (section ^) from perturbation of the mean field 
polaron (thin dashed line) ; 

2. The polaron band is hybridized with the free electron 
band (thin full line) ; 

3. The polaron band is calculated with the TEA approx- 
imation (thick dashed line) ; and 

4. The polaron band is calculated with the HTEA ap- 
proximation where small and large polarons are hy- 
bridized (thick full line) (see next section for details) 

When the quantum lattice fluctuations are small (which 
occurs either at the adiabatic limit a = or when t is 
small), these approximations yield practically the same re- 
sult. When a the best variational form is that which 
gives the lowest energy for the ground-state (that is, the 
largest binding energy). The results of these approxima- 
tions become significantly different when t approaches the 
critical value tp at which the adiabatic first order tran- 
sition occurs. Each of these approximations improves the 
previous one, since the polaron energy becomes lower at 
each step. 

It is clear that approximations (|l|) and (H) , which keep 
the polaron shape rigidly fixed to that at the adiabatic 
limit, are not appropriate to remove the first-order transi- 
tion (see H). The TEA approximation (^) also yields first- 
order transitions, but there are two distinct transitions 
occurring at i = tp{a) < tp and t = tp{a) > tp and the 
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Fig. 5. Binding energy of the quantum polaron versus t at 
a — 0.017 calculated with several approximations as explained 
in the text and magnification (insert). 
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Fig. 6. Variational parameters versus t for the TEA approxi- 
mation of the polaron \(K) (thick lines) and [i(K) (thin lines). 
[i.{K) (thin lines). Wave vector K is zero and a — 0.017 (full 
lines) and a — 0.03 (dashed lines) 



amplitudes of energy discontinuities are weak because the 
polaron shape is determined self consistently. 

Fig.^ shows A and /i values that minimize (^) aX K = 
0. The first of the TEA transitions {i^{a) < tp) occurs be- 
tween a small and a large polaron (see fig|6|) and tp{a) de- 
creases when a increases before that transition disappears 
for a > 0.03. The second TEA transition (ip(a) > tp) per- 
sists for large (a > 0.7) but it is hardly distinguishable on 
the binding energy plot (fig0). The transition occurs be- 
tween a large polaron and a quasi-free electron with an 
extended phonon part (that is, fi{K — 0) tends to 1 when 
t is large). Note also that X{K — 0) may become negative 
in the regime of large t and small a but then the polaron 
binding energy becomes negligible so that it is meaning- 
less to use a polaron picture for a regime that is better 
described as a Fermi liquid. 



4.3 Smoothing the First Order transition: HTEA 

Actually, any first-order transition for the polaron ground- 
state (or the bipolaron) which would be obtained by any 
variational method cannot exist physically. The reason is 
that at the transition point there are two approximate 
wavefunctions with the same variational energy which are 
supposed to approximate the ground-state. It is possible 
to hybridize these two degenerate states to obtain a new 
state with a lower energy. The same arguments hold for 
the exact ground-state, which cannot exhibit any first- 
order transition. 

On the basis of these arguments, we demonstrate nu- 
merically that the two first-order transitions obtained with 
the TEA of polaron can be smoothed using a variational 
form for the wave function !f'^(0) in eq.p3| which hybridizes 
three wave functions, 

tf'^(O) - /3i<(0) + P2l'i{0) + P^'S'ii^)- (43) 

Each of these wave functions has the TEA form ( |4l| ) 
and (^) with parameters Ai,/zi, A2,/i2 and A3, /is re- 
spectively. Hybridizing three wave functions instead of 
two has the advantage of sweeping out simultaneously the 
two successive first order transitions. The variational en- 
ergy ( |3^ ) now depends on 9 parameters Xs,Hs,Ps with 
S e {1,2,3}. Let us note 

Ms.s'iK) =Y^e^^P < ^i{j)\H\^i,{j+p) > 
p 

(44) 

and 

p 

(45) 

where {S,S') G {1,2,3}^. We point out that because of 
the central symmetry of the TEA, the 3x3 matrices M 
and P are real. Then the energy of the ground-state E{K) 
has the following variational form: 

<n-{K)\H\n-{K)>^^^^-^l^l^ (46) 

The extremalization of E{K) eq.^ with respect to /3* 
and (31 yields the set of three equations 

Y.PsMs,s' - E{K){Y,(3sPs.s') = 0, 

s s 

and therefore we have to solve eigenvalue problem M/3 = 
E{K)P(}: that is, E{K) is the lowest eigenvalue of the ma- 
trix P~^/^{K)M{K)P-^/'^. That calculation is very simi- 
lar to the perturbative method of the mean-field described 
previously in the case of the bipolaron but here the lowest 
eigenvalue E{K) has still to be minimized with respect to 
the set of six parameters {Xs,fJ-s)- 

For small t we recover the TEA results (that is, only 
one P is nonnegligible, figj|). Close to the TEA first-order 
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Fig. 7. Binding energy of the quantum polaron versus t cal- 
culated with TEA (thin lines) and the HTEA (thick lines) ( 
a = 0.01). Magnification of the two first order transitions of 
the TEA (inserts). 



transitions the variational groundstate appears as the hy- 
bridization of either a small polaron and a large polaron 
or a large polaron and a quasi-free electron (very large 
polaron). A significant increase of the binding energy of 
polaron results from this hybridization in these crossover 
regions where the first order transitions are smoothed and 
thus removed fig.0. Furthermore our calculations show 
that the energy gain due to hybridization persists for large 
t values. In that regime the fluctuations of the quantum 
lattice are strong enough to hybridize two TEA states, the 
large polaron and the quasi-free electron, whose energies 
differ only slightly. 

A consequence of the hybridization can be also ob- 
served on the shapes of the polaron bands. In the adiabatic 
limit {a — 0), the small polaron is degenerate under ar- 
bitrary lattice translations, which means that the polaron 
band is perfectly flat, as shown on fig.^. In the regime 
where the polaron is metastable for t > tp, the flat po- 
laron band intersects the free electron band so that there 
is a line of wave vectors where the small polaron state 
and the free electron state are degenerate (see fi-g.||) . With 
nonvanishing quantum lattice fluctuations ( a ^ 0), the 
degeneracies are lifted along the intersection line. Approx- 
imation 1^ provides an important correction in the vicin- 
ity of tp where the adiabatic polaron becomes extended. 
Around the degenerate line there is a cross-over region in 
wave vector where the component of the free electron to 
the groundstate varies from almost 1 to almost when K 
goes from to tt and the opposite for the component of 
the small polaron. Thus there is a smooth exchange of the 
quantum state from the large to the small polaron (fig.|lO|). 
This exchange also occurs for the upper band but in re- 
verse order from the small polaron to the large polaron. 

In the band of the TEA, for t(a)^ < t < t{af only 
one first-order transition is observed in K space between 
a small polaron and a large one. For t > t(a)^ two first 
order transition might be observed in K space at different 



W i' 

\\ • I extended 
i V / 1 electron 

V// 




adiabatic | j 
polaron ^ 



Fig. 8. Polaron Bands Es{K) and free electron band El(K) in 
the Kx direction at a = (thin lines) and a = 0.017 calculated 

1) as a perturbation of the adiabatic polaron (thick lines) and 

2) by hybridizing the polaron band and the free electron band 
(thicker lines) at t = 0.0745 < tp (left) and t = 0.085 > tp 
(right). 




Fig. 9. TEA Bands (thin lines) and HTEA Bands (thick fines) 
Es{K) (full lines) and El{K) (dot-dashed) in the Kx direction 
at Q = 0.017 at t = 0.055 < tp{a) (left) and t = 0.06 > tp{a) 
(right) 



K values and they occur between the three kinds of po- 
laron describded previously. They are smoothed with the 
HTEA. 

The effective polaron mass is the inverse curvature 
1/Tp at the bottom of the lowest polaron band at K — 0. 
It can be calculated as a function of t and a and compared 
with the value l/T^p obtained from the Lang-Firsov trans- 
formation fl) . The variation versus t of the ratio Tp/T^p is 
shown fig.nQfor different a. For a small, this ratio Tp/T^p 
is almost one, which confirms that the mean- field approxi- 
mation used to establish formula (^ is valid for both quan- 
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and the extended electron. This sharp variation becomes 
smoother and smoother as a increases. 

Fig.|ll| compares the binding energy of the polaron and 
its tunneUing energy Tp. For small t, the binding energy 
is much larger than Tp, while for large t it becomes much 
smaller. There is a value t — tc{a) where both energies are 
equal. In the vicinity of this region the polaron has max- 
imum mobility while it remains reasonably well-bound 
(compared to this tunnelling energy!). 



5 Variational Calculation of Quantum 
Bipolarons 



Fig. 10. Ratio Tp/Tlf versus t at a = 0.01 (dashed line), 
a — 0.022 (dot-dashed lines), a — 0.03 (full line) calculated as 
a perturbation to the adiabatic limit (thin lines) and using the 
TEA (thick lines). 




Fig. 11. Polaron tunnelling energy Tp (thick lines) and its 
binding energy (thin lines) versus t a.t a = 0.01 (left), a — 0.09 
(right) calculated with the HTEA method. Inserts: Ratio Tp/t 
versus t. 



tum lattice fluctuations and t small. When t increases from 
zero the ratio Tp/Tlf increases from unity, which means 
that eq.^ overestimates the polaron effective mass. We al- 
ready observed this effect for the bipolaron case in the 
formula eq.^. However this effect does not imply sharp 
variations. 

To compare the polaron mass and the bare electronic 
mass, the ratio Tp/t is plotted fig-llf] for different values of 
a. For large t as well as for large a, the polaron effective 
mass reduces to the bare electron mass. In other words 
the electron becomes practically free. 

When a is small, there is a sharp drop in the inverse ef- 
fective polaron mass, which is reminiscent of the first order 
transition at a = between the localized small polaron 



The variational methods ( |34| ) we used for the single po- 
laron can be extended to bipolarons with variational forms 
(SO), (SI) and (QS). For this purpose we write the bipo- 
laron wave function as a Bloch wave: 



3 



(47) 



and we postulate an extended Toyozawa form for the local 
wave function 



|if^(0) > = 



X exp 



> 



(48) 



5.1 TEA quantum Bipolarons 



The simple TEA approximation for the bipolaron consists 
in choosing V'^^, with the form (^ for B=(SO), (|o|) for 
B=(SI) or ( pl| ) for B=(QS) and vi with exponential forms 
which depend on the type of bipolaron as follows: 



..SI 



l.\ + \ly\ 



(0) = -C5i[A*l:i'+""'+A*"^-^l+" 



(0) 



(49) 
(50) 
(51) 



The same arguments used to prove equation (Hn) imply 



(52) 



which determines the parameters Cso, C.qi and Cqs- Us- 
ing the scalar product formula ( p7|) and (gsj), the varia- 
tional energy ( p5[ ) is calculated numerically and minimized 
with respect to both A and /i parameters for each value 
of the wave vector K and for each bipolaron (SO), (SI) or 
(QS). This variational form still has a small number of 
parameters which allows a fast numerical minimization. 
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Fig. 12. TEA energy versus t for bipolaron (SO) (thin full line), 
(SI) (thin dashed line), (QS) (thin dot-dashed line) and energy 
of two single HTEA polarons (full thick line) at i7 = 0.23 and 
a = 0.01. 



Fig. 13. Parameters A (full lines) and (dashed lines) versus 
t for the minima of the variational form ( |35| ) for the TEA 
bipolarons (SO), (SI) and (QS) (from left to right) at (7 = 0.23 
and a = 0.01. Note that (QS) (right figure) has two A variables 
Ai (full line) and A2 (dot-dashed). 



The minimum energy is always found to be at the bot- 
tom of the lowest band at K = 0. The quantum corrections 
to the energies of bipolarons (SO), (SI) and (QS) are com- 
pared with the energy of two polarons far apart. We use 
the HTEA result described in the previous section (fig.p^, 
since we know that it yields the lowest and thus the most 
accurate energy for the quantum polaronic ground-state. 

As for the TEA polaron each TEA bipolarons (SO), 
(SI) or (QS) exhibits a first-order transition when t in- 
creases between a small and a large bipolaron with the 
same symmetry. Actually if one compares the energies of 
all the possible solutions these large bipolarons are found 
never to be the groundstate whatever a is, because a pair 
of single quantum polarons has always less energy. As a 
result, these bipolarons always gain energy by breaking up 
into two polarons (fig, 12) even for large a. 

In the adiabatic limit {a — 0), these TEA calculations 
become identical to the variational calculation which was 
described in (see fig.|^. Comparing the energies of 
these TEA bipolarons (without any hybridization) we con- 
struct a new phase diagram for a non zero with first-order 
transition lines and test how it changes when the quantum 
lattice parameter a increases. 

The approximate calculations of the tunnelling energy 
for the polaron (||) and for the bipolaron (|l^) suggests 
that for U sufficiently different from 1/4 the tunnelling of 
a single polaron with one electron is much easier than that 
of a bipolaron, which contains two electrons and moreover 
involves a bigger lattice distortion. Thus, one should ex- 
pect more generally that the energy gain generated by the 
quantum lattice fluctuations for the single polaron is sys- 
tematically much larger than that for the bipolarons. As 
a result, the domain of parameters where the groundstate 
consists of an unbound pair of large polarons should ex- 
tend at the expense of the domains of the bipolarons when 
a increases. 




Fig. 14. Phase diagram of the TEA bipolarons (SO), (SI), (QS) 
and a pair of unbound polarons for a = (thin dashed line), 
a = 0.001 (thick dashed line), a = 0.0016 (full thick line), and 
a = 0.01 (full thick line?). The case q = is the adiabatic case 
already calculated in ref. 



Indeed fig.|l^ confirms that the first-order transition 
lines which exist in the adiabatic limit shift to lower values 
of t when a increases. As a consequence the domain of 
existence of the (QS) bipolaronic ground-state shrinks to 
zero for a rather small value (approximately 0.002) of a 
and completely disappears for larger values. 

The disappearance of the triple point of the phase dia- 
gram between bipolarons (SO), (SI) and (QS) for relatively 
small values of a seems to rule out our suggestion that 
bipolarons could become very light. We show how to re- 
cover this possibility in the last section by minor changes 
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in the model which may restore this triple point for rela- 
tively large values of a. 

5.2 HTEA quantum Bipolarons 

We carefully examined whether the HTEA calculation of 
bipolarons, which in principle should be more accurate, 
could change this conclusion. Actually, it will not change 
it, and to not confuse or bother the reader, all details of 
our unsuccessful (but useful) numerical investigations are 
not presented. 

As said previously, in principle no first-order transi- 
tions could exist for the ground-state of a pair of elec- 
trons interacting with the lattice. They are removed by 
hybridization of all (or only those which are relevant), de- 
generate bipolaron solutions (SO), (SI) or (QS) both small 
and large, which necessarily generates some energy gain. 

The HTEA calculation for the bipolaron is similar to 
that for the polaron except that it may involve more bipo- 
laronic states. We assume generally that the wave function 
!Z'^(0) (^^ is a normalized combination of n wave func- 
tions (n depends on the number of TEA which hybridize) 

*'^(0)>=^/3sf'|(0)> (53) 
s 

which may have different bipolaronic forms S=(SO), (Sl^;), 
(Sly), (QS) each of which can be small and large, so that 
in principle there are 8 different states. However, we have 
not use simultaneously all these states since there are no 
situations in the parameter space t, U where all their en- 
ergies are simultaneously degenerate but only a relevant 
subset 0. 

Then the energy of the ground-state E{K) has the 
following variational form: 

<n^mH\U^iK)>^ ^^^-^ff;''\ (54) 

1^S.,S' PsPs'Ps,s' 

where 

Ms,s'{K) = < ^sij)\H\^E'{3+p) > 

p 

(55) 

and 

p 

(56) 

The extremalization of (jsj) with respect to (3s is done 
by a diagonahzation of the matrix P^^^'^{K)M{K)P~^/^ 
of size n X n. The variational energy is minimized with 
respect to parameters of eq.(|3|). 

^ Actually using all of them practically does not change the 
result because the irrelevant states hardly hybridize with the 
others. 



In all regions of the phase diagram, for small a the 
HTEA energy corrections for the bipolarons (SO), (SI) or 
(QS) are systematically much smaller than those involved 
by the polarons. The tunnelling energy of bipolarons is 
much smaller than those of the polaron. 

The hybridization cross-overs which are found at each 
smoothed first-order transition of the TEA phase diagram 
remain very narrow and the hybridization energy gain is 
negligible. One needs to have a high bipolaronic degener- 
acy such as the triple point or a relatively large value of 
a (a > 0.05) to observe non negligible crossovers. Even 
in that case the energy gains remain small compared to 
those of an unbound pair of the HTEA polarons. 

If the HTEA bipolarons keep almost the same energy 
as the TEA bipolarons, the phase diagram figjlj is practi- 
cally unchanged. Of course, the first-order transition lines 
which appear in this phase diagram should now be viewed 
as sharp crossover lines. The crossover between the bound 
bipolarons and the unbound pair of polarons has been in- 
vestigated with a general HTEA bipolaron form (including 
the latest) but no significant hybridization has been found 
between these two kinds of states so that we can not draw 
a conclusion about the nature of this transition. 

The triple point is a special point of the phase diagram 
where the bipolarons (SO), (SI) and (QS) are degenerate 
and where we should expect a higher energy gain by hy- 
bridization when a is not too small. Unfortunately, this 
triple point disappears when a increases beyond approxi- 
mately 0.002. When it just disappears the TEA bipolaron 
binding energy referred to two unbound HTEA polaron is 
just zero but then its tunnelling energy T), is maximum 
(but still only 10^^ x t : that is the bipolaron effective 
mass is seven order of magnitude larger than those of the 
bare electron). 

The negative conclusion of this section is that more so- 
phisticated variational calculations does not confirm the 
conclusion of section (^) which was based on the assump- 
tion a small extrapolated to larger a. 

The present study also shows that in the domain of 
small U one may have a quantum bipolaron ground-state 
with a large tunneling energy occurring at very large a > 
0.1. This result is simply obtained with only the TEA of 
the small bipolaron (SO) that is proved to have a nonnegli- 
gible binding energy for both t and U small enough. Nev- 
ertheless, this result is not relevant for such large a, our 
approach based on a perturbative theory of the adiabatic 
limit fails because of too large quantum lattice fluctua- 
tions. 



6 Phonon Dispersion Effect 

We intend to show that highly degenerate point that could 
persist under large quantum lattice fluctuations implies 
very light bipolarons. To achieve that goald, a simple pro- 
cedure consists in changing the model so as to favor the 
bipolaron (QS). If we could make it more robust to quan- 
tum lattice fluctuations it should become very light for 
reasonably large a by hybridization with the other degen- 
erate bipolarons at the triple point. 
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We choose to introduce a phonon dispersion, but this 
might not be the unique way. When an electron is present 
at a given site it will also distort the lattice at the neigh- 
boring sites. If the sign of the dispersion is appropriate, the 
lattice potential at the neighboring sites is lower which fa- 
vors its occupancy by electrons and thus the spatial exten- 
sion of the bipolaron. The bipolaron (QS) which is more 
extended than the bipolaron (SO) should be favored. ^ 

We consider the new Hamiltonian 



Hd^H- c ^ (a+ -I- ai)(a+ + aj) 



(57) 



where Ti. is the Holstein-Hubbard Hamiltonian (|l|) and its 
reduced Hamiltonian corresponds to H (||) that gives 



with 



<ij> 



C = i£(i^)2 
Eo hujo 



(58) 



(59) 
(60) 



When the coupling C is positive the dispersive term 
generates an effective attractive interaction between po- 
larons. This coupling cannot exceed the value 1/2 beyond 
which the low wavevector phonons becomes unstable. 



6.1 Adiabatic Limit 

At the adiabatic limit the equation eq.|l^ becomes 



where D is the matrix 



A,.,: = 1 



A., = 



otherwise 



(62) 



Bipolarons (SO), (SI), (QS)... which were found at the 
anti-integrable limit of the Holstein-Hubbard model at 
t = persist as ground-states in this model with nonzero 
coupling C (see the diagrams figsjl^ and |l^). The do- 
main where bipolaron (QS) is the ground-state enlarges 
when C increases up to its maximum value 1/2. As ex- 
pected the existence of a positive dispersion favors the 
quadrisinglet ground-state. 

The first-order transition between bipolaron (SO) and 
(QS) becomes almost second-order and difficult to dis- 
tinguish numerically since there is no symmetry breaking 

^ Phonon dispersion may induce other important effects in 
the bipolaron structure as shown in ^] for CDW's. 
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Fig. 15. Ground-state phase diagram for Hamiltonian ( |58[ ) 
at C = 0.1 (thick full lines) compared to the initial case C — 
( thin full lines), and approximate diagrams calculated with 
the bipolaron exponential ansatz (thin dashed lines) for same 
couplings C — 0., C — 0.1. 
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Fig. 16. Frequencies versus U of the internal modes of bipo- 
laron (SO) (full lines), (SI) (dashed lines), (QS) (dot-dashed 
lines), for t = 0.04 (left) and t = 0.09 (right). The breathing 
modes are represented by thick lines and the pinning modes 
by thin lines. Vertical lines determine the location of the first 
order transitions. 



between these two bipolaronic states. Then as expected, 
there is a soft internal mode which almost vanishes at the 
transition on both side of the transition which corresponds 
to a breathing mode of the bipolaron with the same sym- 
metry. Simultaneously the Peierls-Nabarro barrier almost 
vanishes. 

This soft mode which does not break the bipolaron 
symmetry is not a pinning mode and does not favor the 
classical mobility of this bipolaron. To that purpose the 
pinning mode which also softens at the first order transi- 
tion between ( QS) and (SI) is the most appropriate ( see 
fig.|l|and refs. ||,i|). 
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Fig. 18. Profiles of electronic density (empty symbol) and absolute value of the displacement (full symbol) versus site in x 
direction for the adiabatic bipolarons (SO), (SI) and (QS) at the triple points: C = 0. t = 0.078 U = 0.233 (squares linked full 
lines), and C = 0.3 t = 0.0904 U = 0.222 (circles linked by dot-dashed lines). 
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Fig. 17. Same as figjl^ but with a phonon dispersion C — 0.3 
calculated exactly in the adiabatic limit (thick full lines) and 
approximated with the exponential ansatz (thin full lines). 



6.2 Quantum Corrections 

Same methods, as those used above for the original Holstein- 
Hubbard model are applied to deal with the quantum lat- 
tice fluctuations of the modified model. The degeneracy 
due to the translation invariance of the model is lifted 
according to standard perturbation theory. One gets a 
tight binding model as in section |^ which yields both 
binding and tunnelling energies of the quantum ground- 
state. Figs.[r^) shows these quantities for a strong coupling 
(C = 0.3). 

The binding energy of the bipolaron refers to two non- 
interacting polarons calculated with the HTEA method, 
which is the most accurate. For a single polaron condition 
(^ becomes 



1 



2(1.-2. *C) 



(63) 



and we choose to write the displacement as 
where v is given by 



D- 



When the coupling C is too large 0.2 < C < 0.5 our ac- 
curacy is limited in practice because of the bipolaron (QS) 
extension, which requires large system sizes we cannot 
afford. This problem occurs especially close to the first- 
order transition between bipolaron (QS) and the extended 
state (see diagram fig.|l^) . However the exponential ansatz 
eq.^ 20 2^ still fits quite well the bipolaron ground-state 
as shown on the diag.^. These variational forms allow us 
to compute quickly the bipolarons even for large C values 
and to determine approximately the groundstate with a 
reasonable accuracy (see diagrams figs.|l^ and p7| ). 

As we already know the flaw of this approximate method 
is that spurious first-order transitions may occur. This sit- 
uation happens nearby the first-order transition between 
(QS) and the extended state as seen on diag.|l^. It is due 
to the exponential ansatz which does not provide a good 
fit of the bipolaron when it becomes more extended. 

However, at the triple point the bipolaron ground-state 
is still localized on very few sites (fig.|l8|) for C = 0.3 and 
the exponential ansatz remains sufficiently accurate. 



|ix| + |iy| 



(64) 



For a large enough phonon coupling C > 0.2, in the 
region we investigate i < 0.1 the HTEA method only 
requires the hybridization between a small polaron and 
a large polaron. The almost second-order transition dis- 
played by the TEA at C = occurs now at a larger tp(a)^. 

The binding energy of the quantum bipolaron is still 
large in that region and one notices the optimal regime 
where both tunnelling and binding energies have the same 
value. 

To obtain the optimal region, a fine tuning of the pa- 
rameters is required because changing them slightly can 
either reduce the binding energy so that the bipolaron be- 
comes fragile against temperature or sharply increase its 
effective mass, killing its quantum mobility. 

Phonon dispersion favors the mobility of the bipolaron 
because it extends the lattice distortion around the bipo- 
laron (see fig.^S]) as well as the electronic wave function. 
Classically, this effect is manifested by internal mode soft- 
ening and by the depression of the Peierls-Nabarro energy 
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barrier (not calculated here see paper I between the 
different bipolarons. As a result, when the lattice is quan- 
tum the hybridization between the different bipolarons is 
increased, which increases the band width and decreases 
the effective mass. 

The HTEA calculatio n (p3|) for the bipolaron confirms 
these properties (see figs.|19D. Condition ( |52| ) becomes 



(1.-2. *C) 



(65) 



and — D~^v where v is still given by eq.(64). 

In the vicinity of the (QS) region (see figs.[l9[) the ef- 
fective mass of the HTEA for the bipolaron is about five 
times larger than the effective mass computed with the 
perturbative method, but the bipolaron mass is still very 
small. The comparison of the binding energy calculated 
with the two methods shows that the variational HTEA 
method is not accurate in the area of the QS region. Indeed 
the perturbative method gives a stronger binding energy 
and thus it is variationally better. This is likely due to the 
fact that when the bipolaron extends too much the TEA 
is not accurate because the bipolaron shape is not well 
approximated by the exponential. 

Fig.|l^ shows for a = 0.017 the effective mass of the 
bipolaron in the optimal regime that ranges not far from 
100 bare electronic mass. We choose as an example the 
realistic optical phonon frequency Tiluq — l.lO^^eV and to 
be in the optimal regime C = 0.3, a = 0.017 U = 0.25 
t = 0.1 the initial parameter of Hamiltonian (^ must be 
g = S.lQ-'^eV Eq = 6eV v = l.5eV t = 0.6eV c = O.SeV. 
The tunnelling energy as well as the bipolaron binding en- 
ergy are about 6.10~'^ey. With such characteristic values 
and a bipolaron concentration not too large, a superfluid 
state could be expected at relatively high temperatures, 
perhaps few hundred degrees K. This estimate neglects 
the bipolaron interactions, but when their concentration 
becomes large these interactions cannot be neglected, es- 
pecially at half filling when there is one polaron per site. 
Close to this close packing regime the bipolaronic struc- 
ture cannot exist anymore for sure. Instead, a magnetic 
spatially ordered polaronic structures could occur. Fur- 
ther studies should investigate the situation with large 
electron densities. 



We realized this situation in a modified Holstein-Hub- 
bard model, which involves both an electron-phonon inter- 
action and a direct repulsive electron-electron interaction. 

The superconducting state of such very light bipo- 
larons occurs for weak concentrations. When the concen- 
tration becomes larger there are strong interactions be- 
tween the bipolarons, which may both break them into 
polarons and organize different structures (for example, 
magnetic) . 

This situation may happen in superconducting cuprates. 
In the undoped regime where the band of electrons is half 
filled, the structure can be viewed as close-packed polarons 
with an antiferromagnetic ordering. This polaron struc- 
ture should persist for low doping till a certain electron 
concentration where the holes are polaron vacancies. For 
a sufficiently large doping the electron concentration may 
become low enough in order that a (first order) transi- 
tion toward a superfiuid of light quantum bipolarons takes 
place. The real phenomenology should be more complex 
because one should expect that the model parameters de- 
pend on the doping and thus that the system does not re- 
main always close from the optimal regime with strongly 
bonded light bipolarons but move around this point. Oth- 
erwise, we suggested in p5[ | that in some appropriate mod- 
els the (QS) bipolaron could have a d-symmetry. We have 
not yet realized an explicit model where such an effect 
occurs, but we hope to. 

The numerical techniques we used (Toyozawa Expo- 
nential Ansatz) and its improvement (HTEA) where the 
hybridization between different states is taken into ac- 
count, turned out to be very efficient to study the bipo- 
laron mass. It should be developed to consider models 
with many electrons. In ||28|| it was proven that at adia- 
batic limit, the ground-state at large electron-phonon cou- 
pling was bipolaronic. Variation of the exponential ansatz 
may provide strong simplifications for these case and a 
qualitative understanding of the many-polaron problem 
first in the adiabatic limit, next with quantum lattice fluc- 
tuations. Finally, the problem of quantization of discrete 
breathers can be approached with similar techniques [p6| . 
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7 Concluding Remarks 

In some circumstances the bipolaron might become un- 
usually light, which allows in principle the formation of 
superconducting states at rather high temperature with 
physically realistic parameters. This effect is due to the 
degeneracy of several bipolaronic states in the adiabatic 
limit for some specific regions of the phase diagram. In 
this situation there are small Peierls Nabarro barriers and 
phonon softening for the different bipolaronic states. Then 
the quantum lattice fiuctuations lift the degeneracy be- 
tween the degenerate states and may yield very light hy- 
bridized bipolarons, which however are well-bound. 
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Abstract. We study the effective mass of the bipolarons and essentially the possibility to get both light 
and strongly bound bipolarons in the Holstein-Hubbard model and some variations in the vicinity of the 
adiabatic limit. Several approaches to investigate the quantum mobility of polarons and bipolarons are 
proposed for this model. First, the quantum fluctuations are treated as perturbations of the mean- field 
(or adiabatic) approximation of the electron-phonon coupling in order to calculate the bipolaron bands, ft 
is found that the bipolaron mass generally remains very large except in the vicinity of the triple point of 
the phase diagram (see j^), where the bipolarons have several degenerate configurations at the adiabatic 
limit (single site (SO), two sites (SI) and quadrisinglet (QS)), while the polarons are much lighter. This 
degeneracy reduces the bipolaron mass significantly. Next we improve this result by variational methods 
(modified Toyozawa Exponential Ansatz or TEA) valid for larger quantum perturbations away from the 
adiabatic limit. We first test this new method for the single polaron. We find that the triple point of the 
phase diagram is washed out by the lattice quantum fluctuations which thus suppress the light bipolarons. 
Further improvements of the method by hybridization of several TEA states do not change this conclusion. 
Next we show that some model variations , for example a phonon dispersion may increase the stability of 
the (QS) bipolaron against the quantum lattice fluctuations. We show that the triple point of the phase 
diagram may be stable to quantum lattice fluctuations and a very sharp mass reduction may occur, leading 
to bipolaron masses of the order of 100 bare electronic mass for realistic parameters. Thus we argue that 
such very light bipolarons could condense as a superconducting state at relatively high temperature when 
their interactions are not too large, that is, their density is small enough. This effect might be relevant for 
understanding the origin of the high Tc superconductivity of doped cuprates far enough from half filling. 

PACS. 71.10.Fd Lattice fermion models (Hubbard model, etc.) - 71.38.-|-i Polarons and electron phonon 
interactions - 74.20.Mn Nonconventional mechanisms (spin fluctuations, polarons and bipolarons, res- 
onating valence bond model, anyon mechanism, marginal Fermi liquid, Luttinger liquid, etc.) - 74.25.Jb 
Electronic structure 

1 Introduction are described by a hard core boson model which could 

become superfluid. Unfortunately, their calculations also 
1.1 Specific Problem for high TC Superconductivity show that when the electron phonon couphng increases 
$H beyond the Migdal instability, the effective mass of these 

- ■ Superconducting materials at temperatures significantly bosons grows exponentially fast and becomes so huge that 
higher than the maximum Tc predicted by MacMiUan ^^eem^s hopdes^s to get supCTConductivity m this model, 
1^ for the standard BCS superconductivity [gj are excep- ^ ^'^^ ^ temperatures, 

tional 0]. Up to now, there is a wide variety of such ma- One of us(SA) conjectured in that the interplay 
terials which are all cuprates built with Cu02 planes and of the electron-phonon coupling with a direct electron- 
with many kinds of interlayer dopants |^ . electron repulsion could reduce significantly the bipolaron 

When the electron-phonon coupling increases too much, effective mass and thus favor high Tc superconductivity, 
it is known for several decades (Migdal |^) that the BCS In the present paper we support this conjecture by cal- 
theory should break down because of lattice instabilities, culating the effective mass of a single bipolaron in a 2d 
These instabilities are associated with the formation of model which involves both an electron-phonon coupling 
polarons and bipolarons. Alexandrov et al. developed and a direct electron-electron repulsion. For this purpose 
later a theory of bipolaronic superconductivity where the we choose to investigate first this effect in the Holstein- 
electrons form strongly bound on-site bipolarons. They Hubbard model because of its simplicity. 
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In the absence of Hubbard repulsion, we confirm that 
the effective mass of the bipolaron is indeed very large 
Ip^ , which is incompatible with a high Tc superconduct- 
ing phase. When the Hubbard term is increased new bipo- 
laronic states become stable. They are 2-site bipolarons 
which consist of two neighboring polarons bound by their 
magnetic interaction in a singlet state and also a bipo- 
laron called a "quadrisinglet" (QS) which consists of the 
combination of four singlets sharing one common site. In 
the parameter region where these bipolaronic states have 
nearly degenerate energies, the effective bipolaron mass 
is sharply reduced. The drastic mass reduction is due to 
resonance between the different bipolarons. Certain vari- 
ations of the model, such as a phonon dispersion, might 
even increase the binding energy of the bipolaron while 
allowing it to keep a very light effective mass for realistic 
parameters. 

We first discuss some early known results about the 
effective masses of polarons and bipolarons. Adiabatic re- 
sults (section |2]) described in are briefly recalled. Next 
we treat the quantum lattice fluctuations as a perturba- 
tion of the adiabatic limit. The main effect is to lift the 
bipolaron degeneracy both due to translational invariance 
and to the possible existence of several kinds of adiabatic 
bipolarons with almost the same energy (section ^ . This 
correction is only valid for a very small quantum lattice 
parameter. 

To extend the field of validity of our calculations we 
next propose to use the Toyozawa variational form. Quan- 
tum polarons and bipolarons are approximated by a self- 
consistent Bloch wave that is exact at the adiabatic limit. 
At this limit it is demonstrated in that there is not a 
great loss of accuracy if the shapes of polaron and bipo- 
laron are exponentials. Thanks to this approximation we 
gain much simplicity for the variational form. We first ap- 
ply this method to polarons (section U) and next to bipo- 
larons (section ^ . For very small quantum lattice fluctua- 
tions the results obtained by perturbation of the adiabatic 
limit are practically recorvered but there are significant 
deviations when the fluctuations become larger. First or- 
der transitions that cannot exist physically are washed out 
by hybridizing several Toyozawa variational forms. Actu- 
ally, the regime where the bipolaron mass is sharply re- 
duced is swallowed up in the domain where the ground- 
state is unbound polarons. We demonstrate in section (||) 
that this undesirable phenomena can be avoided by varia- 
tions ofthe model which increase the stability of the (QS) 
bipolaron (e.g. a phonon dispersion with the appropriate 
sign). 



where j and k represent lattice sites, T is the trans- 
fer integral between nearest neighbor sites < j,k >. The 
electrons are Fermions represented by the standard anti- 
commuting operators Cj'^ and Cj,a at site j with spin 

(T =t or J,. Hi^a- = C^^^Cj^cr and = -t- n^j^. a'^^ and aj 
are standard creation and annihilation boson operators of 
phonons and hujo is the phonon energy of a dispersion-less 
optical phonon branch, g is the on-site electron phonon 
coupling and v the on-site electron-electron repulsion (Hub- 
bard interaction). 

We choose Eq = Sg'^/fiLjQ as energy unit as in 0. 
Defining the position and momentum operators as 



hujQ 



hujo ^ 



with the commutation relation [uj^pj] 
we obtain the dimensionless Hamiltonian 



(2) 

(3) 

(4) 



^^2 



(5) 



<j,k>,cr 



Our reduced dimensionless parameters are 
T „ 1 

t ■ 



Eo = Sg'/hLOo 
T 



7 



1 ?kJOs4 

4 2.9 ' 



(6) 



Despite the primitive nature of our model, it may catch 
important aspects of reality. However, we shall also demon- 
strate at the end of this paper that certain model varia- 
tions could be favorable for bipolaron mass reduction. For 
example, we may introduce a coupling between nearest 
neighbor atoms so that the phonon branch is no longer 
dispersionless. Then the new Hamiltonian is the sum of (|g) 
and the extra energy term —C J2<i j> UiUj where < i,j > 
represents all the pairs of nearest neighbor sites i and 
J. When C > 0, the bipolaron mass reduction is en- 
hanced while it remains strongly bound. Thus we demon- 
strate that relatively minor changes in the model may fa- 
vor (or disfavor) superconductivity at relatively high tem- 
perature. 



1.3 Polaron and Bipolaron Effective Mass 



1.2 The Holstein-Hubbard Model 

Let us first recall our notations for the model we study 
here. Its Hamiltonian is 

+gnj (a+ + a^) + vrij^^Uj^i (1) 



Let us first briefly recall some standard results about the 
effective masses of polarons and bipolarons. The Lang- 
Firsov unitary transformation yields a new Hamilto- 
nian 

3 

<j,k>,cr j 



L. Proville, S. Aubry; Small Bipolarons in the 2-dimensional Holstein-Hubbard Model 



3 



where Hlf = e-"-^^" He''^'^'' with Slf = \ Y.jPj'^j 

After this transformation, the creation operator Cj*"^ 
at site j acts on the vacuum by creating both an electron 
and a lattice distortion 



+ p-»Pj/2| 



> 



(8) 



that is, a polaron. A standard but rough mean-field ap- 
proximation consists in taking the average of the transfer 
integral for unperturbed phonons. We obtain an approxi- 
mate formula for the transfer integral of this polaron 



masses are so huge that it is unreasonable to consider that 
they could exhibit a Bose condensation |]l0| . 

We perform here a numerical calculation of the effec- 
tive mass of the bipolaron (and also the polaron) in order 
to show that in some specific regions of the parameter 
space, when the Hubbard term becomes comparable with 
the electron-phonon binding energy, these effective masses 
can be drastically reduced so that Bose condensations of 
bipolarons become plausible. 



(9) 2 The Mean-field Holstein-Hubbard Model 



When there is a single electron in the system the electron- 
electron interaction does not play any role. Then the ef- 
fective mass of a single polaron is defined as the inverse 
of the second order derivative versus the wave-vector q 
of the polaronic energy E{q), that is, Tlf- The effective 
mass of the polaron is that of the bare electron multiplied 
by exp [l/(8a)]. It becomes exponentially large when a is 
small. 

This approximation tends to become right only when 
the operator corresponding to the transfer integral in eq.^ 
has small fluctuations. This condition is fulfilled when the 
pre-factor t is small. It is also fulfilled when a is large, 
that is, for a weak electron-phonon coupling g compared 
to the phonon energy hiVQ. In the antiadiabatic limit {a 
large), the model becomes a Hubbard model with an on- 
site effective electron-electron interaction U = (C/ — 1/4) 
which is attractive for t7 < 1/4 and repulsive for [/ > 1/4 
and where the transfer integral has been renormalized. 
The negative U model is expected to have superconduct- 
ing phases [|l2[ for non- vanishing band filling. However, we 
treat here the opposite case a small which is close to the 
adiabatic limit. 

Actually when t is small, our numerical results agree 
with formula (^V For larger t, the effective mass of the po- 
laron given by (|9|) becomes larger than the mass we com- 
pute. Note that our result should be more reliable because 
it yields a lower variational energy for a single polaron. 

The effective mass of the bipolaron has been calcu- 
lated by Alexandrov et al in the same limit (t small) 
for strongly bound bipolarons (that is, for U small) and 
considering the kinetic energy term in eq.|^ as a pertur- 
bation. In our dimensionless units they found the transfer 
integral th for a bipolaron 



tb 



At' 



.-l/(2a) 



1-4C/ 



(10) 



If one extrapolates naively this formula for larger U , 
one would find that % becomes infinite. Of course, the 
associated effective mass of the bipolaron cannot vanish, 
but our results nevertheless demonstrate that it is sharply 
depressed not far from the region U « 0.25. Comparison 
of formula (^) and ^ shows that in the region where both 
a and t are small the effective mass, of the bipolaron is 
much larger by many order of magnitude than the polaron 
mass which is itself much larger than that of the bare elec- 
tron. In most physically realistic situations, the bipolaron 



We calculate first the adiabatic bipolarons which are ground- 
state of a mean-field Hamiltonian. They are the exact so- 
lutions in the adiabatic limit when 7 = is zero (that 
is, when the atomic kinetic energy is negligible). These 
spatially localized solutions are degenerate under lattice 
translation. This degeneracy is lifted when the atomic ki- 
netic energy is taken into account. Within a perturbative 
treatment, this explicitly gives bands of extended bipo- 
larons characterized by a wave-vector. The inverse cur- 
vature of the lowest band at zero wave-vector yields the 
bipolaron effective mass. This calculation become exact in 
principle in the limit of small 7. Note that similar meth- 
ods were already developed in ITs] to calculate the effective 
masses of discommensurations in Charge Density Waves. 



2.1 The adiabatic regime 

The Hamiltonian eq.^ can be written as the sum of three 
terms 

H = H,i+Hph + Hf (11) 

where Hei and Hph are decoupled electron and phonon 
Hamiltonians respectively and Hf is a fiuctuation term. 



<i,j>,a- 



i 



(12) 
(13) 

(14) 



fii and Ui are variational parameters which are deter- 
mined by minimizing the ground-state energy of the ef- 
fective Hamiltonian Had — H^i + Hph- It comes out that 
Ui — (ui) and fii — (rii) are the average of the correspond- 
ing operators. The standard mean-field approximation for 
the electron phonon coupling consists in neglecting the 
fiuctuation energy Hf. 

Minimizing the ground-state energy of Hamiltonian 
eqjl^ also yields 

(u,) - -(n,)/2 (15) 
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Then, the ground-state of the mean field Hamiltonian 
Had has the form 



> (16) 



A pair of electrons with the electronic wave function 
ipij is created as well as the corresponding lattice distor- 
tion Ui. The electronic wave function is a singlet state, 
that is, a symmetric function of (i, j): = ipj.i ■ It ful- 
fills an extended nonlinear Schroedinger equation which is 
exactly the same as in the adiabatic case at a = O . 



t 



+ UStj Tp.,J = Eel^puj (17) 



A is the four-dimensional discrete Laplacian and Ui — 
""■1,1 + "-ij with 

n^A^y,H^.J\' (18) 



j 

The square root of the mean square lattice fluctuation 
{{ui — Ui)'^)^^'^ = a = 7^/^ is small of order a. Thus, the 
mean-field approximation obviously becomes exact in the 
adiabatic limit a — > when there are no lattice fluctua- 
tions. 



2.2 Bipolarons from Anti-integrable limit and 
variational approximations 

For an easy understanding the reader should refer to our 
early paper jl] were the adiabatic (or mean field) bipo- 
larons were investigated in detail in the two-dimensional 
model by continuation from the anti-integrable limit {t = 
0)- . 

The main result of m is that we found a quite rich 
phase diagram with first order transition lines in the pa- 
rameter space (U,t). For large t the electrons remain ex- 
tended and do not self localize as bipolarons. For small 
t there are several kinds of structures that compete to 
be the bipolaron ground-state. These bipolarons were de- 
noted (SO), (SI) and (QS). Bipolaron (SO) is mostly local- 
ized at a single site and has square symmetry. Bipolaron 
(SI) consists into a bound pair of polarons in a magnetic 
singlet state localized on two neighboring sites. It breaks 
the square symmetry and is oriented either in the x di- 
rection (Sl)x or the y direction (Sl)j,. The quadrisinglet 
bipolaron (QS) is a combination of four singlet states with 
a common central site and has square symmetry. 

Interesting properties are obtained at a triple point 
corresponding to the intersection of three first-order tran- 
sition lines. At that point, and apart from the transla- 
tional degeneracy, there are four different degenerate bipo- 
larons (SO), (QS), (Sl)^ and (Sl)^. We shall see that the 
quantum lattice perturbations hybridize these degenerate 
states and hence drastically enhance the bandwidth of 
the bipolaron or, equivalently, reduce its effective mass. 
Within a classical picture we already noticed that the en- 
ergy barrier ( Peierls-Nabarro barrier) which has to be 



overcome to move the bipolaron through the lattice was 
drastically reduced. 

We also investigated in [0 some approximations with 
exponential variational forms for the bipolarons that al- 
low analytical calculations. The exact phase diagram cal- 
culated numerically was reproduced with a quite good ac- 
curacy with the following forms 



(19) 



ix-l| + |ial + bx| + |jj,|) 



c 



(20) 



| + |jy±l 



) 



for bipolarons (SO), (SI) and (QS) respectively. A, B, and 
C are normalization factors and the parameters A, Ai and 
A2 are optimized for energy minimization. We shall de- 
velop here a quantum analogous version of these approxi- 
mations to improve our methods. 



3 Quantum Lattice Corrections 

We now treat the mean- field fluctuation Hf = 1/2 X]i("i^ 
Ui){ni — fii) as a perturbation that lifts the translational 
degeneracy of the mean field bipolarons ( [l^ ) , whose wave 
functions are denoted \Q^{j) > where S represents bipo- 
larons (SO), (Sl)a;, (Sl)j^ or (QS). The index j is the site 
where the bipolaron (S) is located. For bipolarons (Sl)a; 
and {Sl)y which occupy two adjacent sites {jx,jy) and 
Ux + l,jy) or {jx,jy) and {jx,jy + 1), respectively, we 
choose by convention j = {jx,jy)- To treat the mean- field 
fluctuation in lowest order, the initial Hamiltonian ( |Tl| ) 
should be projected and diagonalized in the subspace gen- 
erated by all these translated wave functions. We already 
noticed that the bipolaron energies might be degenerate 
or almost degenerate so that we should take into account 
their possible hybridization. The eigenstates should have 
the general form 

l^>=E«^^^l^''0')> (22) 

where as,j are coefficients to be determined by extremal- 
ization of < f2\H\f2 > with the normalization constraint 
< n\n >= 1. Both < n\{2 > and < i7|H|i7 > are 
quadratic functions of asj 

< [2\H\[2 > = E a*s,^M(^s\^),iS',J)as'.J (23) 
<n\Q>= J2 a*s^,P(^sMS',j)as',j (24) 

where the coefficients of matrices P and M are defined as 

P{s,^).iS'.J) = < f2s{^)\ns'U) > (25) 
M(s,,),(S'.j) = < ns{t)\H\ns'{j) > (26) 
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It is important to take into account the fact that the eigen- 
states of the self- consistent mean-field Hamiltonian Had 
are not orthogonal one with each other, since the matrix of 
scalar products is not diagonal. For two normalized wave- 
functions |J7 > and |J7' > with the form ( p^ and with 
electronic wave functions {ipij} and j}, electronic den- 
sities fii = ~2wi and h'^ — —2u[ respectively the scalar 
products eqs.(^6|),(^) can be calculated explicitly for the 
Hamiltonian IM). 




< n\H\f2' > ^ N- < f2\f2' > 

n i,j 
+ i E("" + + '^j.ni'j.n) 

+u J2 v^mV^m - ^ E (^m-^^^.)] (28) 

where ^V'ij' = V'fe,; is the discrete Laplacian on a 4d 
lattice. 

The extremalization equation of < [2\H\fi > with re- 
spect to A = {as,j} with the normalization condition 
A*. P. A = 1, is M.A = E.P.A where E is the Lagrange 
parameter, which is also the eigenenergy. It can be writ- 
ten as an eigenvalue problem for the normalized vector 
B = P^/^A 

P^/^.M.P^/^.B = E.B (29) 

Note that the extensive term < J7|i7' > in the 
second term of eq. (p8|) does not disturb the calculations. 
It yields a constant term in the effective Hamilto- 
nian P~^^^.M.P~^^^ which is nothing but the zero point 
phonon energy of the system with size N (without elec- 
trons). 

Because of the translation invariance of the model, 

M(s,i),(S',j) and_Rs,j),(S'j) O'^ly depends on j - i = n. 
As a result, eq.g9[can be diagonahzed as combinations 
of plane waves with the form Bsj{K) — J2s ^siK)e'^^^ 
with wave vector K which fulfills the eigenequation 

P"^/2(if).M(if).P"^/2(if).B = E^iK).BiK) (30) 
with the Fourier coefficients 

Ps,s'(if)-E^(^^.J%(s'..+«)e'''" (31) 
Ms,s'{K) = E ^'^(5,j),(S',,+„)e'^" (32) 



Then the diagonalization of the 4x4 matrix 
F~^^'^{K).M{K).P~^^'^{K) yields the eigenenergies 
Ei,{K). Figure (|l|) shows an example of calculation of 
these four bands in the vicinity of the triple point. Thus 
when there are four bipolarons that are metastable (e.g. 
in the vicinity of the triple point of the phase diagram) 
one obtains four bipolaron bands E^{K). Within our ap- 
proach the number of bipolaron bands is equal to the 
number of metastable states for the adiabatic bipolaron 
which provides the base about which we expand the eigen- 
states. In other regions of the phase diagram the number 
of metastable bipolarons changes, which induces (unrealis- 
tic) discontinuities for the number of bands. For example, 
when U — only the bipolaron (SO) is metastable, and 
there is only one bipolaron band. 

However the lowest bipolaron band does not exhibit 
very sharp changes despite a small discontinuous varia- 
tion. The reason that the upper bands are not reliable is 
that the energies of these states might be also degener- 
ate with phonon excitations of the bipolaronic states of 
the lower band. The real excited states involve complex 
hybridization between these states. 

Conversely, the bipolaronic states with the lowest en- 
ergies should not hybridize significantly with the higher 
energy states involving phonon excitations. Thus we con- 
sider that the lowest-energy bipolaron band provides a 
reliable description of the bipolaron excitations close to 
its ground-state. We use it to measure the bipolaron ef- 
fective mass, that is, the inverse of the curvature Ti, at 
zero wave vector = 0. Tf, is constant in all directions 
because of the square symmetry (figj^). It can be viewed 
as the effective hopping coefficient for the bipolaron tun- 
nelling through the lattice and can be compared with the 
prediction of given by formula (|o|), which is valid at 
both U and t small. 

Fig.H exhibits the ratio Tf, /if, as a function of the ef- 
fective transfer integral t for U = for several values of 
the quantum parameter a. For U small this ratio goes to 
1 when t goes to 0, which confirms the validity of formula 
( |l0|) predicted by Alexandrov et al. in that regime. We 
also note that beyond this regime when the parameters 
([/, t) are larger than 0, Tf, becomes significantly larger 
than ib, or equivalently the bipolaron effective mass cal- 
culated numerically drops faster than predicted by (^^. 

The insert of fig.^ shows the bipolaron energy gain 
compared to a pair of free electrons |^ at L/ = and the 
fig.y shows the effective transfer integral compared to the 
bare electronic transfer integral, which is negligible at the 
scale of the electronic energy. The bipolaron effective mass 
appears to be much beyond than 10^" electronic masses 
even when its binding energy vanishes. It is clear that 
this regime C/ = is not favorable at all for the Bose 
condensation of such bipolarons that should occur below a 
critical temperature inversely proportional to the effective 
mass of the quasi-particle. 

^ The energy gain compared to a pair of free electrons is 
not an accurate binding energy for the bipolaron. The bipo- 
laron binding energy is precisely measured with respect to an 
unbound polaron state that is defined in the next section. 
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Fig. 2. Ratio Tb/tb versus t of the transfer integral numerically 
calculated and analytically predicted by formula Q {U = 
and a = 0.022 (dot-dashed line), a = 0.017 (dashed line), 
a = 0.01 (full line). 



When the Hubbard term increases for relatively small 
i, fig.§ {t = 0.04) shows that a sharp discontinuity occurs 
when the ground-state bipolaron becomes (SI) U > 0.25. 
There is a sharp increase of the tunnelling energy by five 
orders of magnitude for this bipolaron at a = 0.01. In that 
case the Peierls-Nabarro barrier calculated in the previous 
paper is still very high and consequently there is almost 
no hybridization between (SO) and (SI). The smoothing 
of the discontinuity of the tunnelling energy is thus hardly 
visible. 

When t is larger, the bipolaron (QS) becomes stable for 
U w 0.23 and hybridization between (SO), (SI) and (QS) 
becomes significant. Actually the most important contri- 
bution to the tunnelling energy of the bipolaron comes 
from the hybridization between (QS) and (SI). It is re- 
sponsible for the sharp increase of the tunnelling energy 
or equivalently the sharp drop of the bipolaron effective 
mass. This quantum mobility is favored when (QS) and 
(SI) are degenerate in energy and separate by a weak 
Peierls-Nabarro barrier. Then (QS) may tunnel to one of 
the four neighboring (SI) and the latter tunnels to its 
neighboring (QS) that corresponds to the initial one trans- 
lated by one lattice spacing in the direction of (SI) and so 
on. The bipolaron tunnelling energy could reach 10~'^ the 
bare electronic energy which is not negligible anymore. 

Let us point out that such a high mobility cannot 
be obtained within the approximations used in which 
do not consider the possible degeneracies of several bipo- 
larons. The conclusions of jl^ about the physical impos- 
sibility of bipolaronic superconductivity are irrelevant for 
that situation. 



Fig. 1. Four bipolaron bands E,^{K) versus wave- vector K 
computed close to the triple point where bipolarons (SO), (SI) 
and (QS) are degenerate in energy (a = 0.017, U — 0.232, t — 
0.08). Energies increase from bottom to top. 
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Fig. 3. Ratio Tt/t versus U &t t = 0.04 and a = 0.01; insert: 
Bipolaron energy gain compared to a pair of free electrons 
versus t at U = 0. 



0.27 
U 



Fig. 4. Ratio Tt/t versus f/ at t = 0.04 (top), Tb/t versus U 
at t = 0.08 (bottom) (a = 0.01). 



4 Variational Calculation of Quantum 
Polarons 

In principle the above approach is valid only for very small 
a: that is, when the quantum lattice fluctuations are small. 
However these fluctuations may increase drastically, espe- 
cially close to the first-order transitions when there are 
several degenerate bipolarons that we are especially inter- 
ested in. Thus it is worthwhile to improve our previous 
calculations by a variational approach which should be 
equivalent to the mean field perturbation for small quan- 
tum lattice fluctuations. 



Our purpose is now to develop a quantum version of 
the variational forms used and tested in the adiabatic 
case but which could hold for larger values of a. Our ap- 
proach is a simplified version of those of Toyozawa (see 
[ p^ and [ p^ for a recent application to the polaron in 
ID). We first test this method for the single polaron and 
will extend it in the next section for the bipolarons of the 
Holstein-Hubbard model. 

Because of the invariance of the system under trans- 
lations the wave function of a quantum single polaron is 
written as a Bloch wave: 



\nP{K) >-- 



1 



|<f^(n) > 



(33) 



where yl is a normalization factor and \^^{n) > is ob- 
tained from a unique wave function |!f'^(0) > changing 
all the indices i of its electronic and atomic variables into 
i + n. This transformation is nothing but a shift of the 
wave function from site to n. 



4.1 Toyozawa approximation 

A simple variational approximation proposed by Toyozawa 
is to assume that the local wave function is similar to the 
mean-field polaron: 

E i^k-,C+).e^p (^zY,v,^,.p}j 10 > (34) 

To simplify the spin of the electron is omitted. For each 
wave vector K the variational energy 



(35) 



is a function of the scalar products which does not depend 
on j 

< f^(j)|<f^(i+p) > = exp-^^(w,+p-'y,)^ 



^ (^E^^+p^y 

i 

71 i 



(36) 



(37) 



and has to be extremalized with respect to the 2N param- 
eters corresponding to the electronic wave function {V'J'} 
and the lattice distortion {vi\. This form becomes exact 
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in the adiabatic limit and should improve the previous 
perturbation theory as it is self-consistent. 

A relation between the electronic density and the av- 
erage of the atomic displacement can be easily taken into 
account in this variational form. First, let us recall that the 
true eigenfunctions of iJ are extrema of < > in the 

full space of normalized functions ^ . For a given normal- 
ized eigenfunction f2 of H we can consider the one parame- 
ter family of normalized functions ^{S) — exp (iS.pj)\f2 > 
where the coordinate Uj of the atom j is changed into 
Uj + d. The variational energy of this wave function is 
< 1r{S)\H\1r{S) > which is equal to < 1r{0)\Hs\^{0) > 
where Hs = exp —(ipj6)H enp (ipjS) is simply obtained 
from H by changing Uj into Uj — 5. The variational energy 

< ^{5)\H\^{5) >=< [2\H\[2 > 
- i(5 (2 < n\uj\n > + < {2\nj\n >) + ^6^ (38) 

should be extremal for 6 = 0, which implies 

< n\uj\n >= < n\nj\n > . (39) 

This result is nothing but an extension to the non-adiabatic 
case of the standard relation between the average atomic 
positions and the electronic densities. 

If we now consider an extremum of < f2\H\f2 > for 
I n > normalized in the variational space defined by eq. ( |3^ ) 
and (|3^) this space is no longer invariant under the unitary 
operator e'xjp{iSpj}, but it still remains globally invariant 
under operator exp{i(5 pj } which performs a uniform 
displacement by 5 of all the atoms. We apply the same 
argument as above that is, study < n{S)\H\n{S) > where 
f2{5) = exp{i5J2jPj}^ is extremal for S = 0. This con- 
dition yields J2j < ^Wjl^^ >= -1/2 < n\nj\0 >. 
For the variational extrema with the Toyozawa form ( ^3|) 
and m, we find < '?'(0I E„ W"|>f (m) >= (E„««) < 
'P{l)\<Pjrn) >, which readily implies J^j < > = 

X^n'*^"- ^'^^ polaron, that is, for a system with only 
one electron the extremum of the Toyozawa form ( p3| ) and 
( ^ ) necessarily fulfills 

E^" = -^ (40) 

n 



4.2 Toyozawa Exponential Approximation: TEA 

Minimizing the variational form (^) for the whole set of 
2N — 1 parameters {ipi} and {vi} with condition ( pO| ) is a 
complex numerical task which moreover will become even 
more complex when extended to the bipolaron problem. 
However, we can expect that the behavior of the varia- 
tional parameters {(j^n} and will not be far from ex- 
ponential at infinity. Thus assuming simple exponentials 
for {0„} and Jw„} should not be a bad approximation as 
proposed in at the adiabatic limit. Taking into ac- 
count the normalization and condition (Eoh we postulate 



that the electronic wave function and the atomic modula- 
tion have the form: 

V^f = AAI^-I+I*"! = (1 + A2)/(1- A^) (41) 

v[ = -Bf,\^^\+^^y\ B-' = 2{^ + ^l)/{l-^,) (42) 

for each wave vector K there are only two variational 
parameters X{K) and ii{K) instead of 27V for the original 
Toyozawa ansatz which allows much simpler calculations 
although we still need a numerical minimization of (p5|). 
To that aim we use a simplex method p^ , which is the 
most efficient algorithm we tested because it avoids any 
precision problem due to the numerical computation of 
the derivatives. 

The Toyozawa Exponential Ansatz (TEA) turns out 
to be almost as good as the full ansatz when the polarons 
are small, since in that case the exponential approximates 
quite well its shape. When the size of the polaron becomes 
larger, the TEA (as well as the original Toyozawa ansatz) 
yields a first order transition. This first order transition 
is well-known to exist at the adiabatic limit a = where 
the ground state of a single electron undergoes a first order 
transition from a small polaron to a free electron [p^p^ 
att^tpfn 0.07486. 

We define the binding energy of the quantum polaron 
as the difference between the energy of the extended elec- 
tron at zero wave vector K = and that of the bottom of 
the polaron band. Fig|| shows the variation of the bind- 
ing energy versus t for the quantum polaron calculated in 
several different approximations including the assumption 
that: 

1. The polaron band is calculated as for the bipolaron 
bands (section ^) from perturbation of the mean field 
polaron (thin dashed line) ; 

2. The polaron band is hybridized with the free electron 
band (thin full line) ; 

3. The polaron band is calculated with the TEA approx- 
imation (thick dashed line) ; and 

4. The polaron band is calculated with the HTEA ap- 
proximation where small and large polarons are hy- 
bridized (thick full line) (see next section for details) 

When the quantum lattice fluctuations are small (which 
occurs either at the adiabatic limit a = or when t is 
small), these approximations yield practically the same re- 
sult. When a the best variational form is that which 
gives the lowest energy for the ground-state (that is, the 
largest binding energy). The results of these approxima- 
tions become significantly different when t approaches the 
critical value tp at which the adiabatic first order tran- 
sition occurs. Each of these approximations improves the 
previous one, since the polaron energy becomes lower at 
each step. 

It is clear that approximations (|l|) and (H) , which keep 
the polaron shape rigidly fixed to that at the adiabatic 
limit, are not appropriate to remove the first-order transi- 
tion (see H). The TEA approximation (^) also yields first- 
order transitions, but there are two distinct transitions 
occurring at i = tp{a) < tp and t = tp{a) > tp and the 
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Fig. 5. Binding energy of the quantum polaron versus t at 
a — 0.017 calculated with several approximations as explained 
in the text and magnification (insert). 
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Fig. 6. Variational parameters versus t for the TEA approxi- 
mation of the polaron \(K) (thick lines) and [i(K) (thin lines). 
[i.{K) (thin lines). Wave vector K is zero and a — 0.017 (full 
lines) and a — 0.03 (dashed lines) 



amplitudes of energy discontinuities are weak because the 
polaron shape is determined self consistently. 

Fig.^ shows A and /i values that minimize (^) aX K = 
0. The first of the TEA transitions {i^{a) < tp) occurs be- 
tween a small and a large polaron (see fig|6|) and tp{a) de- 
creases when a increases before that transition disappears 
for a > 0.03. The second TEA transition (ip(a) > tp) per- 
sists for large (a > 0.7) but it is hardly distinguishable on 
the binding energy plot (fig0). The transition occurs be- 
tween a large polaron and a quasi-free electron with an 
extended phonon part (that is, fi{K — 0) tends to 1 when 
t is large). Note also that X{K — 0) may become negative 
in the regime of large t and small a but then the polaron 
binding energy becomes negligible so that it is meaning- 
less to use a polaron picture for a regime that is better 
described as a Fermi liquid. 



4.3 Smoothing the First Order transition: HTEA 

Actually, any first-order transition for the polaron ground- 
state (or the bipolaron) which would be obtained by any 
variational method cannot exist physically. The reason is 
that at the transition point there are two approximate 
wavefunctions with the same variational energy which are 
supposed to approximate the ground-state. It is possible 
to hybridize these two degenerate states to obtain a new 
state with a lower energy. The same arguments hold for 
the exact ground-state, which cannot exhibit any first- 
order transition. 

On the basis of these arguments, we demonstrate nu- 
merically that the two first-order transitions obtained with 
the TEA of polaron can be smoothed using a variational 
form for the wave function !f'^(0) in eq.p3| which hybridizes 
three wave functions, 

tf'^(O) - /3i<(0) + P2l'i{0) + P^'S'ii^)- (43) 

Each of these wave functions has the TEA form ( |4l| ) 
and (^) with parameters Ai,/zi, A2,/i2 and A3, /is re- 
spectively. Hybridizing three wave functions instead of 
two has the advantage of sweeping out simultaneously the 
two successive first order transitions. The variational en- 
ergy ( |3^ ) now depends on 9 parameters Xs,Hs,Ps with 
S e {1,2,3}. Let us note 

Ms.s'iK) =Y^e^^P < ^i{j)\H\^i,{j+p) > 
p 

(44) 

and 

p 

(45) 

where {S,S') G {1,2,3}^. We point out that because of 
the central symmetry of the TEA, the 3x3 matrices M 
and P are real. Then the energy of the ground-state E{K) 
has the following variational form: 

<n-{K)\H\n-{K)>^^^^-^l^l^ (46) 

The extremalization of E{K) eq.^ with respect to /3* 
and (31 yields the set of three equations 

Y.PsMs,s' - E{K){Y,(3sPs.s') = 0, 

s s 

and therefore we have to solve eigenvalue problem M/3 = 
E{K)P(}: that is, E{K) is the lowest eigenvalue of the ma- 
trix P~^/^{K)M{K)P-^/'^. That calculation is very simi- 
lar to the perturbative method of the mean-field described 
previously in the case of the bipolaron but here the lowest 
eigenvalue E{K) has still to be minimized with respect to 
the set of six parameters {Xs,fJ-s)- 

For small t we recover the TEA results (that is, only 
one P is nonnegligible, figj|). Close to the TEA first-order 
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Fig. 7. Binding energy of the quantum polaron versus t cal- 
culated with TEA (thin lines) and the HTEA (thick lines) ( 
a = 0.01). Magnification of the two first order transitions of 
the TEA (inserts). 



transitions the variational ground-state appears as the hy- 
bridization of either a small polaron and a large polaron 
or a large polaron and a quasi-free electron (very large 
polaron). A significant increase of the binding energy of 
polaron results from this hybridization in these crossover 
regions where the first order transitions are smoothed and 
thus removed fig.0. Furthermore our calculations show 
that the energy gain due to hybridization persists for large 
t values. In that regime the fluctuations of the quantum 
lattice are strong enough to hybridize two TEA states, the 
large polaron and the quasi-free electron, whose energies 
differ only slightly. 

A consequence of the hybridization can be also ob- 
served on the shapes of the polaron bands. In the adiabatic 
limit {a — 0), the small polaron is degenerate under ar- 
bitrary lattice translations, which means that the polaron 
band is perfectly flat, as shown on fig.^. In the regime 
where the polaron is metastable for t > tp, the flat po- 
laron band intersects the free electron band so that there 
is a line of wave vectors where the small polaron state 
and the free electron state are degenerate (see fig.||) . With 
nonvanishing quantum lattice fluctuations ( a ^ 0), the 
degeneracies are lifted along the intersection line. Approx- 
imation 1^ provides an important correction in the vicin- 
ity of tp where the adiabatic polaron becomes extended. 
Around the degenerate line there is a cross-over region in 
wave vector where the component of the free electron to 
the ground-state varies from almost 1 to almost when 
K goes from to tt and the opposite for the component of 
the small polaron. Thus there is a smooth exchange of the 
quantum state from the large to the small polaron (fig.^ . 
This exchange also occurs for the upper band but in re- 
verse order from the small polaron to the large polaron. 

In the band of the TEA, for t^{a) < t < t^{a) only 
one first-order transition is observed in K space between 
a small polaron and a large one. For t > t^(a) two first 
order transition might be observed in K space at different 



W i' 

\\ • I extended 
i V / 1 electron 

V// 




adiabatic | j 
polaron ^ 



Fig. 8. Polaron Bands Es{K) and free electron band El(K) in 
the Kx direction at a = (thin lines) and a = 0.017 calculated 

1) as a perturbation of the adiabatic polaron (thick lines) and 

2) by hybridizing the polaron band and the free electron band 
(thicker lines) at t = 0.0745 < tp (left) and t = 0.085 > tp 
(right). 




Fig. 9. TEA Bands (thin lines) and HTEA Bands (thick fines) 
Es{K) (full lines) and El{K) (dot-dashed) in the Kx direction 
at Q = 0.017 at t = 0.055 < (left) and t = 0.06 > 4(q) 

(right) 



K values and they occur between the three kinds of po- 
laron describded previously. They are smoothed with the 
HTEA. 

The effective polaron mass is the inverse curvature 
1/Tp at the bottom of the lowest polaron band at K — 0. 
It can be calculated as a function of t and a and compared 
with the value l/T^p obtained from the Lang-Firsov trans- 
formation fl) . The variation versus t of the ratio Tp/T^p is 
shown fig.nQfor different a. For a small, this ratio Tp/T^p 
is almost one, which confirms that the mean- field approxi- 
mation used to establish formula (^ is valid for both quan- 
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and the extended electron. This sharp variation becomes 
smoother and smoother as a increases. 

Fig.|ll| compares the binding energy of the polaron and 
its tunneUing energy Tp. For small the binding energy 
is much larger than Tp, while for large t it becomes much 
smaller. There is a value t — tc{a) where both energies are 
equal. In the vicinity of this region the polaron has max- 
imum mobility while it remains reasonably well-bound 
(compared to this tunnelling energy!). 



5 Variational Calculation of Quantum 
Bipolarons 



Fig. 10. Ratio Tp/Tlf versus t at a = 0.01 (dashed line), 
a — 0.022 (dot-dashed lines), a — 0.03 (full line) calculated as 
a perturbation to the adiabatic limit (thin lines) and using the 
TEA (thick lines). 



The variational methods (|3J) we used for the single po- 
laron can be extended to bipolarons with variational forms 
(SO), (SI) and (QS). For this purpose we write the bipo- 
laron wave function as a Bloch wave: 




Fig. 11. Polaron tunnelling energy Tp (thick lines) and its 
binding energy (thin lines) versus t a.t a = 0.01 (left), a — 0.09 
(right) calculated with the HTEA method. Inserts: Ratio Tp/t 
versus t. 



tum lattice fluctuations and t small. When t increases from 
zero the ratio Tp/Tlf increases from unity, which means 
that eq.^ overestimates the polaron effective mass. We al- 
ready observed this effect for the bipolaron case in the 
formula eq.^. However this effect does not imply sharp 
variations. 

To compare the polaron mass and the bare electronic 
mass, the ratio Tp/t is plotted fig.|ll] for different values of 
a. For large t as well as for large a, the polaron effective 
mass reduces to the bare electron mass. In other words 
the electron becomes practically free. 

When a is small, there is a sharp drop in the inverse ef- 
fective polaron mass, which is reminiscent of the first order 
transition at a = between the localized small polaron 



and we postulate an extended Toyozawa form for the local 
wave function 




5.1 TEA quantum Bipolarons 

The simple TEA approximation for the bipolaron consists 
in choosing V'^j. with the form (^ for B=(SO), (|o|) for 
B=(S1) or ( pl| ) for B=(QS) and vi with exponential forms 
which depend on the type of bipolaron as follows: 



.f°(0) 



\l.-l\ + \ly 



(0) 



(49) 
(50) 
(51) 



The same arguments used to prove equation (^) imply 



(52) 



which determines the parameters Cso, Csi and Cqs- Us- 
ing the scalar product formula ( p7|) and (p8|), the varia- 
tional energy ( |35|) is calculated numerically and minimized 
with respect to both A and ^ parameters for each value 
of the wave vector K and for each bipolaron (SO), (SI) or 
(QS) (see fig.[l^) . This variational form still has a small 
number of parameters which allows a fast numerical min- 
imization. 
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Fig. 12. TEA energy versus t for bipolaron (SO) (thin full line), 
(SI) (thin dashed line), (QS) (thin dot-dashed line) and energy 
of two single HTEA polarons (full thick line) at i7 = 0.23 and 
a = 0.01. 



Fig. 13. Parameters A (full lines) and (dashed lines) versus 
t for the minima of the variational form ( |35| ) for the TEA 
bipolarons (SO), (SI) and (QS) (from left to right) at U = 0.23, 
a — 0.01 and the wave vector K = 0. Note that (QS) (right 
figure) has two A variables Ai (full line) and A2 (dot-dashed). 



The minimum energy is always found to be at the bot- 
tom of the lowest band at K = 0. The quantum corrections 
to the energies of bipolarons (SO), (SI) and (QS) are com- 
pared with the energy of two polarons far apart. We use 
the HTEA result described in the previous section (fig.p^, 
since we know that it yields the lowest and thus the most 
accurate energy for the quantum polaronic ground-state. 

As for the TEA polaron each TEA bipolarons (SO), 
(SI) or (QS) exhibits a first-order transition when t in- 
creases between a small and a large bipolaron with the 
same symmetry. Actually if one compares the energies of 
all the possible solutions these large bipolarons are found 
never to be the ground-state whatever a is, because a pair 
of single quantum polarons has always less energy. As a 
result, these bipolarons always gain energy by breaking up 
into two polarons (fig, 12) even for large a. 

In the adiabatic limit (a = 0), these TEA calcula- 
tions become identical to the variational calculation which 
was described in 0| (see flg.p^. Comparing the energies 
of these TEA bipolarons (without any hybridization) we 
construct a new phase diagram for a non zero with first- 
order transition lines and test how it changes when the 
quantum lattice parameter a increases. 

The approximate calculations of the tunnelling energy 
for the polaron (||) and for the bipolaron (|l^) suggests 
that for U sufficiently different from 1/4 the tunnelling of 
a single polaron with one electron is much easier than that 
of a bipolaron, which contains two electrons and moreover 
involves a bigger lattice distortion. Thus, one should ex- 
pect more generally that the energy gain generated by the 
quantum lattice fluctuations for the single polaron is sys- 
tematically much larger than that for the bipolarons. As a 
result, the domain of parameters where the ground-state 
consists of an unbound pair of large polarons should ex- 
tend at the expense of the domains of the bipolarons when 
a increases. 




Fig. 14. Phase diagram of the TEA bipolarons (SO), (SI), (QS) 
and a pair of unbound polarons for a = (thin dashed line), 
a = 0.001 (thick dashed line), a = 0.0016 (fuU thin line), and 
a = 0.01 (full thick line). The case a = is the adiabatic case 
already calculated in ref. |l[. 



Indeed fig.|l^ confirms that the first-order transition 
lines which exist in the adiabatic limit shift to lower values 
of t when a increases. As a consequence the domain of 
existence of the (QS) bipolaronic ground-state shrinks to 
zero for a rather small value (approximately 0.002) of a 
and completely disappears for larger values. 

The disappearance of the triple point of the phase dia- 
gram between bipolarons (SO), (SI) and (QS) for relatively 
small values of a seems to rule out our suggestion that 
bipolarons could become very light. We show how to re- 
cover this possibility in the last section by minor changes 
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in the model which may restore this triple point for rela- 
tively large values of a. 

5.2 HTEA quantum Bipolarons 

We carefully examined whether the HTEA calculation of 
bipolarons, which in principle should be more accurate, 
could change this conclusion. Actually, it will not change 
it, and to not confuse or bother the reader, all details of 
our unsuccessful (but useful) numerical investigations are 
not presented. 

As said previously, in principle no first-order transi- 
tions could exist for the ground-state of a pair of elec- 
trons interacting with the lattice. They are removed by 
hybridization of all (or only those which are relevant), de- 
generate bipolaron solutions (SO), (SI) or (QS) both small 
and large, which necessarily generates some energy gain. 

The HTEA calculation for the bipolaron is similar to 
that for the polaron except that it may involve more bipo- 
laronic states. We assume generally that the wave function 
!Z'^(0) (^^ is a normalized combination of n wave func- 
tions (n depends on the number of TEA which hybridize) 

*'^(0)>=^/3sf'|(0)> (53) 
s 

which may have different bipolaronic forms S=(SO), (Sl^;), 
(Sly), (QS) each of which can be small and large, so that 
in principle there are 8 different states. However, we have 
not use simultaneously all these states since there are no 
situations in the parameter space t, U where all their en- 
ergies are simultaneously degenerate but only a relevant 
subset 0. 

Then the energy of the ground-state E{K) has the 
following variational form: 

<n^mH\U^iK)>^ ^^^-^ff;''\ (54) 

1^S.,S' PsPs'Ps,s' 

where 

Ms,s'{K) = < ^sij)\H\^E'{3+p) > 

p 

(55) 

and 

p 

(56) 

The extremalization of (jsj) with respect to (3s is done 
by a diagonahzation of the matrix P^^^'^{K)M{K)P~^/^ 
of size n X n. The variational energy is minimized with 
respect to parameters of eq.(|3|). 

^ Actually using all of them practically does not change the 
result because the irrelevant states hardly hybridize with the 
others. 



In all regions of the phase diagram, for small a the 
HTEA energy corrections for the bipolarons (SO), (SI) or 
(QS) are systematically much smaller than those involved 
by the polarons. The tunnelling energy of bipolarons is 
much smaller than those of the polaron. 

The hybridization cross-overs which are found at each 
smoothed first-order transition of the TEA phase diagram 
remain very narrow and the hybridization energy gain is 
negligible. One needs to have a high bipolaronic degener- 
acy such as the triple point or a relatively large value of 
a (a > 0.05) to observe non negligible crossovers. Even 
in that case the energy gains remain small compared to 
those of an unbound pair of the HTEA polarons. 

If the HTEA bipolarons keep almost the same energy 
as the TEA bipolarons, the phase diagram figjlj is practi- 
cally unchanged. Of course, the first-order transition lines 
which appear in this phase diagram should now be viewed 
as sharp crossover lines. The crossover between the bound 
bipolarons and the unbound pair of polarons has been in- 
vestigated with a general HTEA bipolaron form (including 
the latest) but no significant hybridization has been found 
between these two kinds of states so that we can not draw 
a conclusion about the nature of this transition. 

The triple point is a special point of the phase diagram 
where the bipolarons (SO), (SI) and (QS) are degenerate 
and where we should expect a higher energy gain by hy- 
bridization when a is not too small. Unfortunately, this 
triple point disappears when a increases beyond approxi- 
mately 0.002. When it just disappears the TEA bipolaron 
binding energy referred to two unbound HTEA polaron is 
just zero but then its tunnelling energy T), is maximum 
(but still only 10^^ x t : that is the bipolaron effective 
mass is seven order of magnitude larger than those of the 
bare electron). 

The negative conclusion of this section is that more so- 
phisticated variational calculations does not confirm the 
conclusion of section (^) which was based on the assump- 
tion a small extrapolated to larger a. 

The present study also shows that in the domain of 
small U one may have a quantum bipolaron ground-state 
with a large tunneling energy occurring at very large a > 
0.1. This result is simply obtained with only the TEA of 
the small bipolaron (SO) that is proved to have a nonnegli- 
gible binding energy for both t and U small enough. Nev- 
ertheless, this result is not relevant for such large a, our 
approach based on a perturbative theory of the adiabatic 
limit fails because of too large quantum lattice fluctua- 
tions. 



6 Phonon Dispersion Effect 

We intend to show that highly degenerate point that could 
persist under large quantum lattice fluctuations implies 
very light bipolarons. To achieve that goald, a simple pro- 
cedure consists in changing the model so as to favor the 
bipolaron (QS). If we could make it more robust to quan- 
tum lattice fluctuations it should become very light for 
reasonably large a by hybridization with the other degen- 
erate bipolarons at the triple point. 
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We choose to introduce a phonon dispersion, but this 
might not be the unique way. When an electron is present 
at a given site it will also distort the lattice at the neigh- 
boring sites. If the sign of the dispersion is appropriate, the 
lattice potential at the neighboring sites is lower which fa- 
vors its occupancy by electrons and thus the spatial exten- 
sion of the bipolaron. The bipolaron (QS) which is more 
extended than the bipolaron (SO) should be favored. ^ 

We consider the new Hamiltonian 



Hd^H- c ^ (a+ -I- ai)(a+ + aj) 



(57) 



where Ti. is the Holstein-Hubbard Hamiltonian (|l|) and its 
reduced Hamiltonian corresponds to H (||) that gives 



with 



<ij> 



C = i£(i^)2 
Eo hujo 



(58) 



(59) 
(60) 



When the coupling C is positive the dispersive term 
generates an effective attractive interaction between po- 
larons. This coupling cannot exceed the value 1/2 beyond 
which the low wavevector phonons becomes unstable. 



6.1 Adiabatic Limit 

At the adiabatic limit the equation eq.|l^ becomes 



where D is the matrix 



A,.,: = 1 



A., = 



otherwise 



(62) 



Bipolarons (SO), (SI), (QS)... which were found at the 
anti-integrable limit of the Holstein-Hubbard model at 
t = persist as ground-states in this model with nonzero 
coupling C ||l| (see the diagrams figs, 15 and [l7|) . The do- 
main where bipolaron (QS) is the ground-state enlarges 
when C increases up to its maximum value 1/2. As ex- 
pected the existence of a positive dispersion favors the 
quadrisinglet ground-state. 

The first-order transition between bipolaron (SO) and 
(QS) becomes almost second-order and difficult to dis- 
tinguish numerically since there is no symmetry breaking 

^ Phonon dispersion may induce other important effects in 
the bipolaron structure as shown in jl^] for CDW's. 
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Fig. 15. ground-state phase diagram for Hamiltonian (jS^) at 
C = 0.1 (thick full lines) compared to the initial case C = 
( thin full lines), and approximate diagrams calculated with 
the bipolaron exponential ansatz (thin dashed lines) for same 
couplings C = 0., C — 0.1. 
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Fig. 16. Frequencies versus U of the internal modes of bipo- 
laron (SO) (full lines), (SI) (dashed lines), (QS) (dot-dashed 
lines), for t = 0.04 (left) and t = 0.09 (right). The breathing 
modes are represented by thick lines and the pinning modes 
by thin lines. Vertical lines determine the location of the first 
order transitions. 



between these two bipolaronic states. Then as expected, 
there is a soft internal mode which almost vanishes at the 
transition on both side of the transition which corresponds 
to a breathing mode of the bipolaron with the same sym- 
metry. Simultaneously the Peierls-Nabarro barrier almost 
vanishes. 

This soft mode which does not break the bipolaron 
symmetry is not a pinning mode and does not favor the 
classical mobility of this bipolaron. To that purpose the 
pinning mode which also softens at the first order transi- 
tion between ( QS) and (SI) is the most appropriate ( see 
fig.|l|and refs. ig,|l). 
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Fig. 18. Profiles of electronic density (empty symbol) and absolute value of the displacement (full symbol) versus site in x 
direction for the adiabatic bipolarons (SO), (SI) and (QS) at the triple points: C = 0. t = 0.078 U = 0.233 (squares linked by 
full lines), and C = 0.3 t = 0.0904 U = 0.222 (circles hnked by dot-dashed lines). 
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Fig. 17. Same as figjl^ but with a phonon dispersion C — 0.3 
calculated exactly in the adiabatic limit (thick full lines) and 
approximated with the exponential ansatz (thin full lines). 



6.2 Quantum Corrections 

Same methods, as those used above for the original Holstein- 
Hubbard model are applied to deal with the quantum lat- 
tice fluctuations of the modified model. The degeneracy 
due to the translation invariance of the model is lifted 
according to standard perturbation theory. One gets a 
tight binding model as in section |^ which yields both 
binding and tunnelling energies of the quantum ground- 
state. Figs.[r^) shows these quantities for a strong coupling 
(C = 0.3). 

The binding energy of the bipolaron refers to two non- 
interacting polarons calculated with the HTEA method, 
which is the most accurate. For a single polaron condition 
(^ becomes 



1 



2(1.-2. *C) 



(63) 



and we choose to write the displacement as 
where v is given by 



D- 



When the coupling C is too large 0.2 < C < 0.5 our ac- 
curacy is limited in practice because of the bipolaron (QS) 
extension, which requires large system sizes we cannot 
afford. This problem occurs especially close to the first- 
order transition between bipolaron (QS) and the extended 
state (see diagram fig.|l^) . However the exponential ansatz 
eq.^ 20 2^ still fits quite well the bipolaron ground-state 
as shown on the diag.^. These variational forms allow us 
to compute quickly the bipolarons even for large C values 
and to determine approximately the ground-state with a 
reasonable accuracy (see diagrams figs.|l^ and p7| ). 

As we already know the flaw of this approximate method 
is that spurious first-order transitions may occur. This sit- 
uation happens nearby the first-order transition between 
(QS) and the extended state as seen on diag.|l^. It is due 
to the exponential ansatz which does not provide a good 
fit of the bipolaron when it becomes more extended. 

However, at the triple point the bipolaron ground-state 
is still localized on very few sites (fig.|l8|) for C = 0.3 and 
the exponential ansatz remains sufficiently accurate. 



|ix| + |iy| 



(64) 



For a large enough phonon coupling C > 0.2, in the 
region we investigate i < 0.1 the HTEA method only 
requires the hybridization between a small polaron and 
a large polaron. The almost second-order transition dis- 
played by the TEA at C = occurs now at a larger tp{a). 

The binding energy of the quantum bipolaron is still 
large in that region and one notices the optimal regime 
where both tunnelling and binding energies have the same 
value. 

To obtain the optimal region, a fine tuning of the pa- 
rameters is required because changing them slightly can 
either reduce the binding energy so that the bipolaron be- 
comes fragile against temperature or sharply increase its 
effective mass, killing its quantum mobility. 

Phonon dispersion favors the mobility of the bipolaron 
because it extends the lattice distortion around the bipo- 
laron (see fig.^S]) as well as the electronic wave function. 
Classically, this effect is manifested by internal mode soft- 
ening and by the depression of the Peierls-Nabarro energy 
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barrier (not calculated here see paper I between the 
different bipolarons. As a result, when the lattice is quan- 
tum the hybridization between the different bipolarons is 
increased, which increases the band width and decreases 
the effective mass. 

The HTEA calculatio n (p3|) for the bipolaron confirms 
these properties (see figs.|19D. Condition ( |52| ) becomes 



(1.-2. *C) 



(65) 



and — D~^v where v is still given by eq.(64). 

In the vicinity of the (QS) region (see figs.[l9[) the ef- 
fective mass of the HTEA for the bipolaron is about five 
times larger than the effective mass computed with the 
perturbative method, but the bipolaron mass is still very 
small. The comparison of the binding energy calculated 
with the two methods shows that the variational HTEA 
method is not accurate in the area of the QS region. Indeed 
the perturbative method gives a stronger binding energy 
and thus it is variationally better. This is likely due to the 
fact that when the bipolaron extends too much the TEA 
is not accurate because the bipolaron shape is not well 
approximated by the exponential. 

Fig.|l^ shows for a = 0.017 the effective mass of the 
bipolaron in the optimal regime that ranges not far from 
100 bare electronic mass. We choose as an example the 
realistic optical phonon frequency Tiluq — l.lO^^eV and to 
be in the optimal regime C = 0.3, a = 0.017 U = 0.25 
t = 0.1 the initial parameter of Hamiltonian (^ must be 
g = S.lQ-'^eV Eq = 6eV v = l.5eV t = 0.6eV c = O.SeV. 
The tunnelling energy as well as the bipolaron binding en- 
ergy are about 6.10~'^ey. With such characteristic values 
and a bipolaron concentration not too large, a superfluid 
state could be expected at relatively high temperatures, 
perhaps few hundred degrees K. This estimate neglects 
the bipolaron interactions, but when their concentration 
becomes large these interactions cannot be neglected, es- 
pecially at half filling when there is one polaron per site. 
Close to this close packing regime the bipolaronic struc- 
ture cannot exist anymore for sure. Instead, a magnetic 
spatially ordered polaronic structures could occur. Fur- 
ther studies should investigate the situation with large 
electron densities. 



7 Concluding Remarks 

In some circumstances the bipolaron might become un- 
usually light, which allows in principle the formation of 
superconducting states at rather high temperature with 
physically realistic parameters. This effect is due to the 
degeneracy of several bipolaronic states in the adiabatic 
limit for some specific regions of the phase diagram. In 
this situation there are small Peierls Nabarro barriers and 
phonon softening for the different bipolaronic states. Then 
the quantum lattice fluctuations lift the degeneracy be- 
tween the degenerate states and may yield very light hy- 
bridized bipolarons, which however are well-bound. 



We realized this situation in a modified Holstein-Hub- 
bard model, which involves both an electron-phonon inter- 
action and a direct repulsive electron-electron interaction. 

The superconducting state of such very light bipo- 
larons occurs for weak concentrations. When the concen- 
tration becomes larger there are strong interactions be- 
tween the bipolarons, which may both break them into 
polarons and organize different structures (for example, 
magnetic) . 

This situation may happen in superconducting cuprates. 
In the undoped regime where the band of electrons is half 
filled, the structure can be viewed as close-packed polarons 
with an antiferromagnetic ordering. This polaron struc- 
ture should persist for low doping till a certain electron 
concentration where the holes are polaron vacancies. For 
a sufficiently large doping the electron concentration may 
become low enough in order that a (first order) transi- 
tion toward a superfluid of light quantum bipolarons takes 
place. The real phenomenology should be more complex 
because one should expect that the model parameters de- 
pend on the doping and thus that the system does not re- 
main always close from the optimal regime with strongly 
bonded light bipolarons but move around this point. Oth- 
erwise, we suggested in that in some appropriate mod- 
els the (QS) bipolaron could have a d-symmetry. We have 
not yet realized an explicit model where such an effect 
occurs, but we hope to. 

The numerical techniques we used (Toyozawa Expo- 
nential Ansatz) and its improvement (HTEA) where the 
hybridization between different states is taken into ac- 
count, turned out to be very efficient to study the bipo- 
laron mass. It should be developed to consider models 
with many electrons. In | |2l[ | it was proven that at adia- 
batic limit, the ground-state at large electron-phonon cou- 
pling was bipolaronic. Variation of the exponential ansatz 
may provide strong simplifications for these case and a 
qualitative understanding of the many-polaron problem 
first in the adiabatic limit, next with quantum lattice fluc- 
tuations. Finally, the problem of quantization of discrete 
breathers can be approached with similar techniques [EH . 
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